Accepted by The Astrophysical Journal 

Preprint typeset using IATgX style emulatcapj v. 04/03/99 



HOW EXACTLY DID THE UNIVERSE BECOME NEUTRAL? 



o 

Q 

o 



> 

00 

(N 

On 
On 

Oh! 

6 

in 
c3 



Sara Seager 1 , Dimitar D. Sasselov, 

Astronomy Department, Harvard University, 
60 Garden Street, Cambridge, MA, 02138 
AND 

Douglas Scott 

Department of Physics & Astronomy, 
University of British Columbia, Vancouver, BC, V6T 1Z1 
Accepted by The Astrophysical Journal 

ABSTRACT 

We present a refined treatment of H, He I, and He II recombination in the early Universe. The 
difference from previous calculations is that we use multi-level atoms and evolve the population of each 
level with redshift by including all bound-bound and bound-free transitions. In this framework we follow 
several hundred atomic energy levels for H, He I, and He II combined. 

The main improvements of this method over previous recombination calculations are: (1) allowing 
excited atomic level populations to depart from an equilibrium distribution; (2) replacing the total 
recombination coefficient with recombination to and photoionization from each level calculated directly 
at each redshift step; and (3) correct treatment of the He I atom, including the triplet and singlet states. 

We find that the ionization fraction x c = n c /riH is approximately 10% smaller at redshifts <J 800 than 
in previous calculations, due to the non-equilibrium of the excited states of H, which is caused by the 
strong but cool radiation field at those redshifts. In addition we find that He I recombination is delayed 
compared with previous calculations, and occurs only just before H recombination. These changes in 
turn can affect the predicted power spectrum of microwave anisotropies at the few percent level. 

Other improvements such as including molecular and ionic species of H, including complete heating 
and cooling terms for the evolution of the matter temperature, including collisional rates, and including 
feedback of the secondary spectral distortions on the radiation field, produce negligible change to x e . 
The lower x c at low z found in this work affects the abundances of H molecular and ionic species by 
10-25%. However this difference is probably not larger than other uncertainties in the reaction rates. 

Subject headings: cosmology: theory — atomic processes — early universe — cosmic microwave 
background 



1. INTRODUCTION 

The photons that Penzias & Wilson (1965) detected 
coming from all directions with a temperature of about 
3 K have traveled freely since their last Thomson scatter- 
ing, when the Universe became cool enough for the ions 
and electrons to form neutral atoms. During this recom- 
bination epoch, the opacity dropped precipitously, matter 
and radiation were decoupled, and the anisotropies of the 
Cosmic Microwave Background radiation (CMB) were es- 
sentially frozen in. These anisotropies of the CMB have 
now been detected on a range of scales (e.g. White, Scott, 
& Silk 1994; Smoot & Scott 1998) and developments in 
the field have been very rapid. Recently, two post-COBE 
missions, the Microwave Anisotropy Probe (MAP) and the 
Planck satellite, were approved with the major science goal 
of determining the shape of the power spectrum of anis- 
otropics with experimental precision at a level similar to 
current theoretical predictions. 

Detailed understanding of the recombination process is 
crucial for modeling the power spectrum of CMB aniso- 
tropies. Since the seminal work of the late 1960s (Pee- 
bles 1968; Zel'dovich et al. 1968), several refinements have 
been introduced by, for example, Matsuda, Sato, & Takeda 
(1971), Zabotin & Nasel'skii (1982), Lyubarsky & Sun- 



yaev (1983), Jones & Wyse (1985), Krolik (1990), and 
others, but in fact little has changed (a fairly compre- 
hensive overview of earlier work on recombination is to 
be found in Section IIIC of Hu et al. 1995, hereafter 
HSSW). More recently, refinements have been made inde- 
pendently in the radiative transfer to calculate secondary 
spectral distortions (Dell' Antonio & Rybicki 1993; Rybicki 
& Deir Antonio 1994), and in the chemistry (Stancil, Lepp, 
& Dalgarno 1996b). These improvements may have notice- 
able effects (at the 1% level) on the calculated shapes of 
the power spectrum of anisotropies. Given the potential 
to measure important cosmological parameters with MAP 
and Planck (e.g. Jungman et al. 1996; Bond, Efstathiou, & 
Tegmark 1997; Zaldarriaga, Spergel, & Seljak 1997; Eisen- 
stein, Hu, & Tegmark 1998; Bond & Efstathiou 1998) it is 
of great interest to make a complete and detailed calcula- 
tion of the process of recombination. Our view is that this 
is in principle such a simple process that it should be so 
well understood that it could never affect the parameter 
estimation endeavor. 

Our motivation is to carry out a 'modern' calculation 
of the cosmic recombination process. All of the physics 
is well understood, and so it is surprising that cosmolo- 
gists have not moved much beyond the solution of a sin- 
gle ODE, as introduced in the late 1960s. With today's 
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computing power, there is no need to make the sweep- 
ing approximations that were implemented 30 years ago. 
We therefore attempt to calculate to as great an extent as 
possible the full recombination problem. The other differ- 
ence compared with three decades ago is that we are now 
concerned with high precision calculations, because of the 
imminent prospect of high fidelity data. 

ft is our intention to present here a coupled treatment 
of the non-equilibrium radiative transfer and the detailed 
chemistry. The present investigation was motivated by in- 
dications (HSSW) that multi-level non-equilibrium effects 
in H and He, as well as in some molecular species, may 
have measurable effects on the power spectrum of CMB 
anisotropics by affecting the low-z and high-z tails of the 
visibility function e~ T dT / dz (where r is the optical depth) . 

To that effect this paper presents a study of the recom- 
bination era by evolving neutral and ionized species of H 
and He, and molecular species of H simultaneously with 
the matter temperature. We believe our work represents 
the most accurate picture to date of how exactly the Uni- 
verse as a whole became neutral. 

2. BASIC THEORY 
2.1. The Cosmological Picture 

We will assume that we live in a homogeneous, expand- 
ing Universe within the context of the canonical Hot Big 
Bang paradigm. The general picture is that at some suf- 
ficiently early time the Universe can be regarded as an 
expanding plasma of hydrogen plus some helium, with 
around 10 9 photons per baryon and perhaps some non- 
baryonic matter. As it expanded and cooled there came a 
time when the protons were able to keep hold of the elec- 
trons and the Universe became neutral. This is the period 
of cosmic recombination. 

In cosmology it is standard to use redshift z as a time 
coordinate, so that high redshift represents earlier times. 
Explicitly, a(t) — 1/(1 + z) is the scale factor of the Uni- 
verse, normalized to be unity today, and with the rela- 
tionship between scale factor and time depending on the 
particular cosmological model. It is sufficient to consider 
only hydrogen and helium recombination, since the other 
elements exist in minute amounts. The relevant range of 
redshift is then <^ 10,000, during which the densities are 
typical of those astrophysicists deal with every day, and 
the temperatures are low enough that there are no rela- 
tivistic effects. 

2.2. The Radiation Field 

In describing the radiation field and its interaction with 
matter, we must use a specific form of the radiative trans- 
fer equation. It should describe how radiation is absorbed, 
emitted and scattered as it passes through matter in a 
medium which is homogeneous, isotropic, infinite, and ex- 
panding. The basic time-dependent form of the equation 
of transfer is 

1 dl(f, n, t) dl(f,n,h',t) 

c dt + di = 
j(r, h, v, t) - n{r, h, v, t)I{r, h, v, t). (1) 

Here the symbols are: I(f, n, v, t) is the specific intensity 
of radiation at position r, traveling in direction h (the 



unit direction vector), with frequency v at a time t (in 
units of ergs s _1 cm _2 Hz _1 sr _1 ); I is the path length along 
the ray (and is a coordinate- independent path length); j 
is the emissivity, which is calculated by summing prod- 
ucts of upper excitation state populations and transition 
probabilities over all relevant processes that can release a 
photon at frequency v, including electron scattering; and 
k is the absorption coefficient, which is the product of an 
atomic absorption cross section and the number density of 
absorbers summed over all states that can interact with 
photons of frequency v. 

In the homogeneous and isotropic medium of the early 
Universe, we can integrate equation (Q) over all solid an- 
gles ui (i.e. integrating over the unit direction vector n): 



4?r dJ(r, u, t) 



V--F(r>,f) 



c dt 

j(f,n,v,t) — K(r,h,v,t)I(f,n,v,t) du>. (2) 

Here J(v, t) is the mean intensity, the zeroth order mo- 
ment of the specific intensity over all angles (in units of 
ergs s~ 1 cm _2 Hz _1 ), F is the flux of radiation, which is 
the net rate of radiant energy flow across an arbitrarily 
oriented surface per unit time and frequency interval, and 
c is the speed of light. If the radiation field is isotropic 
there is a ray-by-ray cancellation in the net energy trans- 
port across a surface and the net flux is zero. Also, because 
of the isotropy of the radiation field and the medium being 
static, we can drop the dependence upon angle of j and k 
in equation (|J). With the definition of J, § Idw — 
this simplifies equation (0) to: 



1 dJ(v, t) 
c dT~ 



t) - n(v, t)J(u, t). 



(3) 



The above equation is for a static medium. An isotropi- 
cally expanding medium would reduce the number density 
of photons due to the expanding volume, and reduce their 
frequencies due to redshifting. The term due to the density 
change will be simply a 3^j|y factor, while the redshifting 
term will involve the frequency derivative of J(v, t) and 
hence a v^t^ factor. 

Then the equation for the evolution of the radiation field 
as affected by the expansion and the sources and sinks of 
radiation becomes 



dJ{v,t) dJ{v,t) 



uH(t) 



8J{y,t) 



dt dt dv 

= -ZH{t)J{v, t) + c t) - k{u, t)J{v, t)] , (4) 

where H(t) = a/a. 

This equation is in its most general form and difficult 
to solve; fortunately we can make two significant simpli- 
fications, because the primary spectral distortions 2 are of 
negligible intensity (DelP Antonio & Rybicki 1993) and the 
quasi-static solution for spectral line profiles is valid (Ry- 
bicki & Dell' Antonio 1994). The first simplification is that 
for the purposes of this paper (in which we do not study 
secondary spectral distortions), we set J(v.t) = B(v,t), 
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the Planck function, which is observed to approximate 
J{v) to at least 1 part in 10 4 (Fixsen et al. 1996). Thus we 
eliminate explicit frequency integration from the simulta- 
neous integration of all equ ation s (§2.6). The validity of 
this assumption is shown in §3.5 where we follow the dom- 
inant secondary distortions of H Ly a and H 2-photon by 
including their feedback on the recombination process, and 
find that the secondary spectral distortions from the other 
Lyman lines and He are not strong enough to feed back 
on the recombination process. 

The second simplification is in the treatment of the evo- 
lution of the resonance lines (Lya, etc.) which must still 
be treated explicitly - because of cosmological redshifting 
they cause J(u,t) ^ B{v,t) in the lines. These we call 
the primary distortions. We use esca pe pr obab ility meth- 
ods for moving (expanding) media (§ |2.3.3 and 3.1). This 
simplification is not an approximation, but is an exact 
treatment - a simple solution to the multi-level radiative 
transfer problem afforded by the physics of the expanding 
early Universe. 

Not only does using B{v,t) with the escape probabil- 
ity method instead of J(i/, t) simplify the calculation and 
reduce computing time enormously, but the effects from 
following the actual radiation field will be small compared 
to the main improvements of our recombination calcula- 
tion which are the level-by-level treatment of H, He I and 
He II, calculating recombination directly, and the correct 
treatment of He I triplet and singlet states. 

2.3. The Rate Equations 

The species we evolve in the expanding Universe are 
H I, H II, He I , He II, He III, e~, H~, H 2 , and H 2 + . 
The chemistry of the early Universe involves the reactions 
of association and dissociation among these species, facili- 
tated by interactions with the radiation field, J(i>, t). The 
rate equations for an atomic system with N energy levels 
can be described as 



a(t) 



^d( ni (t)a(t) 3 ) 
alt 



n e (t)n c (t)P ci - m{t)Pk 



+Ef=i n 3 {t)P 3l - ni {t)P l3 , (5) 



where the Pij are the rate coefficients between bound lev- 



ij — j-hj ~r n c Cij and 
n c Cj C , where R refers to radiative rates and C 



els i and j, and the Pj C are the rate coefficients between 
bound levels and the continuum c: Pa — R 

P lc = R_ 

to collisional rates. Here the ns are physical (as opposed to 
comoving) number densities: n, refers to number density 
of the ith excited atomic state, n c to the number density 
of electrons, and n c to the number density of a continuum 
particle such as a proton, He II or He III. a(t) is the cos- 
mological scale factor. The rate equations for molecules 
take a slightly different form because their formation and 
destruction depend s on the rate coefficients for the reac- 
tions discussed in §2.3.4, and molecular bound states are 
not included. 

2.3.1. Photoionization and Photorecombination 

By calculating photorecombination rates R C i directly to 
each level for multi-level H, He I, and He II atoms, we 
avoid the problem of finding an accurate recombination 
coefficient, the choice of which has a large effect on the 
power spectra (HSSW). 



Photoionization rates are calculated by integrals of the 
incident radiation field J(v, t) and the bound-free cross 
section o~i C (y). The photoionization rate in s _1 is 



Ri, 



An 



J{v, t)dv. 



(6) 



Here i refers to the ith excited state and c refers to the con- 
tinuum, vq is the threshold frequency for ionization from 
the ith excited state. The radiation field J(v, t) depends 
on frequency v and time t. With m as the number density 
of the ith excited state, the number of photoionizations 
per unit volume per unit time (hereafter photoionization 
rate) is n,iRi c . 

By using the principle of detailed balance in the case 
of local thermodynamic equilibrium (LTE), the radiative 
recombination rate can be calculated from the photoion- 
ization rate. Then, as described below, the photorecom- 
bination rate can be generalized to the non-LTE case by 
scaling the LTE populations with the actual populations 
and substituting the actual radiation field for the LTE 
radiation field. In LTE the radiation field J(y, t) is the 
Planck function B(y,t). B(v,t) is a function of time t 
during recombination because Tr = 2.728(1 + z(t)) K. We 
will call the LTE temperature T (T = T R = T M at early 
times), where Tr is the radiation temperature, and Tm the 
matter temperature. To emphasize the Planck function's 
dependence on temperature, we will use B(v,T) where T 
is a function of time. 

By detailed balance in LTE we have 



[n e n c R ci ] 



LTE 



WiiRi 



iLTE 



Radiative recombination includes spontaneous and stim- 
ulated recombination, so we must rewrite the above equa- 
tion as 



r D iLTE r D spon-iLTE r D stim~| LTE tn\ 

[n c n c R ci \ = [n c n c R^ \ + [n c n c R ci \ (7) 

fr u iLTE r nstimlL™^ . r nstiml LTE 

= [\niRic\ ~ [mRic \ ) + [n l R lc \ 
Using the definition of i?j C in equation (Q) , 

[n c n c R cl ] LTE = 47m^ E H T) (l - dv 
,ltb/ a ^ v l B ( v ,T) e - h ^l^ T dv. (8) 



47m- 



The first term on the right hand side is the spontaneous 
recombination rate and the second term on the right hand 
side is the stimulated recombination rate. Here hp is 
Planck's constant and /cb is Boltzmann's constant. The 
factor (1 — e ~' lpI/ / fcBT ) is the correction for stimulated re- 
combination (see Mihalas 1978, §4-3 for a derivation of 
this factor). Stimulated recombination can be treated as 
either negative ionization or as positive recombination; the 
physics is the same (see Seager & Sasselov, in prepara- 
tion, for some subtleties). With the LTE expression for 
recombination (equation (||)), it is easy to generalize to the 
non-LTE case, considering spontaneous and stimulated re- 
combination separately. Because the matter temperature 
Tm and the radiation temperature Tr differ at low z, it is 
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important to understand how recombination depends on 
each of these separately. 

Spontaneous recombination involves a free electron but 
its calculation requires no knowledge of the local radiation 
field, because the photon energy is derived from the elec- 
tron's kinetic energy. In other words, whether or not LTE 
is valid, the LTE spontaneous recombination rate holds 
per ion, as long as the velocity distribution is Maxwellian. 
The local Planck function (as representing the Maxwell 
distribution) depends on Tm, because the Maxwell distri- 
bution describes a collisional process. Furthermore, since 
the Maxwellian distribution depends on Tm, so does the 
spontaneous rate. To get the non-LTE rate, we only have 
to rescale the LTE ion density to the actual ion density: 



n c n c 



(n e n c ) 



LTE 



hpv 



D 



n c n c 



LTE 



(i/,T M )(l-e- hp "/ fcB7 ^)di/ (9) 



a ic (u) 2h P is 3 c _ hpu/kBTM 
hpis c 2 



dv. 
( 10 ) 

To generalize the stimulated recombination rate from 
the LTE rate to the non-LTE rate, we rescale the LTE ion 
density to the actual ion density, and replace the LTE radi- 
ation field by the actual radiation field J(v, t) because that 
is what is 'stimulating' the recombination. The correction 
for stimulated recombination depends on Tm, because the 
recombination process is collisional; the term always re- 
mains in the LTE form because equation (^) was derived 
from detailed balance. So we have 



n c n c R^ m = 4tt 



n c n c 



(n e n c ) 



LTE 



r^ TE x 



^lj(v,t)e- h ^' k ^dv. 
hpv 



(11) 



Therefore, the total non-LTE recombination rate 

R^r) is 



spoil 



Tl c Tl c R cl 



4tt 



hpv 



n e n c 
2hpv z 



LTE 



+ J{y, t) 



The LTE population ratios (rii / n c n c ) LTE depend only on 
Tm through the Saha relation: 



n c n c 



LTE 



2-Km c kp,T\ 



M 



3/2 



2g c 



(13) 



Here m c is the electron mass, the atomic parameter g is the 
degeneracy of the energy level, and Ei is the ionization en- 
ergy of level i. In the recombination calculation presented 
in this paper we use the Planck function B(y^ Tr) instead 
of the radiation field J(v, T) as described earlier. 

In the early Universe Case B recombination is used. 
This excludes recombinations to the ground state and con- 
siders the Lyman lines to be optically thick. An implied 
assumption necessary to compute the photoionization rate 
is that the excited states (n > 2) are in equilibrium with 



the radiation. Our approach is more general than Case 
B, because we don't consider the Lyman lines to be op- 
tically thick and don't assume equilibrium among the ex- 
cited states. For more d etails on the validity of Case B 
recombination see §3.2.5, To get the total recombination 
coefficient, we sum over captures to all excited levels above 
the ground state. 

To summarize, the form of the total photoionization rate 



is 



N 

riiRic 

i>l 



N c 
^n;47T / 
i>l Jv 



Q"»c(^) 

hpv 



B(u,T R )du, (14) 



and the total recombination rate is 



N N 



i>l 
4-7T 



LTE 



hpv 



B{u,T R ) 



-h P v/kvT Mdv _ ( 15 ) 



2.3.2. Comparison With the 'Standard' Recombination 
Calculation of Hydrogen 

The 'standard' recombination calculation refers to the 
calculation widely used today and first derived by Pee- 
bles (1968, 1993) and Zel'dovich and collaborators (1968, 
1983), updated with the most recent para met ers and re- 
combination coefficient (HSSW). See also §|3.l| . 

For a 300-level H atom in our new recombination calcu- 
lation, the expressions (Q) and ([l5]) include 300 integrals 
at each redshift step. The standard recombination calcu- 
lation does not go through this time-consuming task but 
avoids it entirely by using a pre-calculated recombination 
coefficient that is a single expression dependent on Tm 
only. The recombination coefficient to each excited state 
i is defined by 

a t (T M ) = RT n - (16) 

Here a is a function of Tm, because spontaneous recombi- 
nation is a collisional process, as described previously. The 
total Case B recombination coefficient (ae) is obtained 
from 

N 

Q!b(Tm) = y^Qj»(T M ). (17) 
i>i 

We will refer to q;b(Tm) as the 'pre-calculated recombina- 
tion coefficient' because the recombination to each atomic 
level i and the summation over i are pre-calculated for 
LTE conditions. See Hummer (1994) for an example of 
how these recombination coefficients are calculated. Some 
more elaborate derivations of a(Tvi) = /(Tm,?i) have also 
been tried (e.g. Boschan & Biltzinger 1998). 

The standard recombination calculation uses a pho- 
toionization coefficient /3b (Tm) which is derived from de- 
tailed balance using the recombination rate: 



(riiPi) 



LTE 



(n c n p ai) 



LTE 



(18) 



To get the non-LTE rate, one uses the actual populations 

Hi, 

LTE 



nij3i(T M ) = rii 



II; 



cti(T M ), 



(19) 
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or with the Saha relation (equation (|l3|)), 



ni(3i{Ti 



Ml 



( 2Ttm Q k B T] 



hi 



3/2 



9i 



(20) 

Constants and variables are as described before. To get 
the total photoionization rate, (3i is summed over all ex- 
cited levels. Because the 'standard' calculation avoids use 
of all levels i explicitly, the m are assumed to be in equi- 
librium with the radiation, and thus can be related to the 
first excited state number density ni s by the Boltzmann 
relation, 

ni ^ n2s 9^ e ^-E t )/ kB T M ^ 
92s 

With this relation, the total photoionization rate is 



(21) 



JV 
i>l 



2'Km c k B T\ 



M 



3/2 



(22) 

In this expression for the total photoionization rate, the 
excited states are populated according to a Boltzmann dis- 
tribution. Tm is used instead of Tr, because the Saha and 
Boltzmann equilibrium used in the derivation are colli- 
sional processes. The expression says nothing about the 
excited levels being in equilibrium with the continuum, be- 
cause the actual values of n e , ri\ and «2s are used, and the 
rii are proportional to n,2 S . 

To summarize, the standard calculation uses a single ex- 
pression for each of the total recombination rate and the 
total photoionization rate that is dependent on Tm only. 
The total photoionization rate is 



^ riiRic = ti2 S Pb(Tm) 



i>l 



n 2s a B {T M )e- E ^l k * T » (2irm c k B T M ) 3/2 (23) 



and the total recombination rate is 



N 

E 

i>l 



n c n p R c 



i p a B (Ti 



MJ 



(24) 



Comparing the right hand side of equation (^3|) to our 
level- by-level total photoionization rate (equation (113)), 
the main improvement in our method over the standard 
one is clear: we use the actual excited level populations rii, 
assuming no equilibrium distribution among them. In this 
way we can test the validity of the equilibrium assump- 
tion. Far less important is that the standard recombi- 
nation treatment cannot distinguish between Tr and Tm, 
even though photoionization and stimulated recombina- 
tion are functions of Tr while spontaneous recombination 
is a function of Tm, as shown in equations ( |i"4| ) and ([l5|). 
The non-equilibrium of excited states is important at the 
10% level in the residual ionization fraction for z < 800, 
while using Tr in photoionization and photoexcitation is 
only important at the few percent level for z < 300 (for 
typical cosmological models). Note that although the pre- 
calculated recombination coefficient includes spontaneous 
recombination only, stimulated recombination (as a func- 
tion of Tm) is still included as negative photoionization via 
detailed balance (see equation (R)). 



2.3.3. Photoexcitation 

In the expanding Universe, redshifting of the photons 
must be taken into account (see equation (Q)). Line pho- 
tons emitted at one position may be redshifted out of in- 
teraction frequency (redshifted more than the width of 
the line) by the time they reach another position in the 
flow. We use the Sobolev escape probability to account for 
this, a method which was first used for the expanding Uni- 
verse by Dell' Antonio and Rybicki (1993). The Sobolev es- 
cape probability (Sobolev 1946), also sometimes called the 
large-velocity gradient approximation, is not an approxi- 
mation but is an exact, simple solution to the multi-level 
radiative transfer in the case of a large velocity gradient. 
It is this solution which allows the explicit inclusion of the 
line distortions to the radiation field - without it our de- 
tailed approach to the recombination problem would be 
intractable. We will call the net bound-bound rate for 
each line transition ARji , where j is the upper level and i 
is the lower level: 

ARji = pij {n.j [Aji + BjiB(uij,t)] - niBijB{vij,t)} . 

(25) 

Here the terms Aji, Bji, B^ are the Einstein coefficients; 
the escape probability pij is the probability that photons 
associated with this transition will 'escape' without be- 
ing further scattered or absorbed. If = 1 the photons 
produced in the line transition escape to infinity - they 
contribute no distortion to the radiation field. If = 
no photons escape to infinity; all of them get reabsorbed 
and the line is optically thick. This is the case of primary 
distortions to the radiation field, and the Planck function 
cannot be used for the line radiation. In general pij <C 1 
for the Lyman lines and pij = 1 for all other line transi- 
tions. With this method we have described the redshifting 
of photons through the resonance lines and found a simple 
solution to the radiative transfer problem for all bound- 
bound transitions. The rest of this section is devoted to 
deriving pij . 

For the case of no cosmological redshifting, the radiative 
rates per cm 3 for transitions between excited states of an 
atom are 



where 



Hi Rij — rii Bij J 

and n j R, j ^ — n j A. j i ~\~ n j B jiJ, 



J = J(v, t)4>(v)dv, 
Jo 



(26) 
(27) 



(28) 



and 0(z^) is the line profile function with its area normal- 
ized by 



(v)dv = 1. 



(29) 



The line profile (j>(v) is taken to be a Voigt function that 
includes natural and Doppler broadening. In principle, 
equation ( p8| ) is the correct approach. In practise we take 
4>{v) as a delta function, and use J{v, t) instead of J. 
The smooth radiation field is essentially constant over the 
width of the line and so the line shape is not important; 
we get the same results using J or J(v,t). 
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The Sobolev escape probability considers the distance 
over which the expansion of the medium induces a veloc- 
ity difference equal to the thermal velocity (for the case of 
a Doppler width): L = Uth/I^'li where t>th is the thermal 
velocity width and v' the velocity gradient. The theory 
is valid when this distance L is much smaller than typical 
scales of macroscopic variation of other quantities. 

We follow Rybicki (1984) in the derivation of the 
Sobolev escape probability. The general definition of es- 
cape probability is given by the exponential extinction law, 



Pij = exp[-r(j/y)], 



(30) 



where i/y is the frequency for a given line transition, and 
T(vij) is the monochromatic optical depth forward along 
a ray from a given point to the boundary of the medium. 
Here 7"(i/y) is defined by 



dr{vij) = -k(j>(i>ij)dl, 



(31) 



where k is the integrated line absorption coefficient, so 
that the monochromatic absorption coefficient or opacity 
is k = k(f>{yij), and I is the distance along the ray from the 
emission point (1 = 0). Rewriting the optical depth for a 
line profile function (which has units of inverse frequency) 
of the dimcnsionlcss frequency variable x — (y — Vij)/A, 
with A the width of the line in Doppler units, and Vij the 
central line frequency, we have 



dr(i/y) 



4>(x)dl. 



(32) 



Here k is the absorption coefficient (defined in equa- 
tion (jl])), with k just dividing out the line profile function, 
and 

k = -^—(riiBij - UjBji). (33) 
Using the Einstein relations giBij — gjBji and Aji — 



(2hv 3 /c 2 )Bji the absorption coefficient can also be written 



as 



i gj 
n 4 — - iij 



(34) 



Note that distances along a ray I correspond to shifts 
in frequency x. This is due to the Doppler effect induced 
by the velocity gradient and is the essence of the Sobolev 
escape probability approach. For example, the Lya pho- 
tons cannot be reabsorbed in the Ly a line if they redshift 
out of the frequency interaction range. This case will hap- 
pen at some frequency x, or at some distance I from the 
photon emission point, where because of the expansion the 
photons have redshifted out of the frequency interaction 
range. For an expanding medium with a constant velocity 
gradient v' = dv/dl, the escape probability along a ray is 
then 



Pij = exp 



k 
'A 



4>{x-l/L)dl 







= exp 


-75 / 4>(x)dx 




J — oo 




(35) 



The velocity field has in effect introduced an intrinsic es- 
cape mechanism for photons; beyond the interaction limit 
with a given atomic transition, the photons can no longer 



be absorbed by the material, even if it is of infinite ex- 
tent, but escape freely to infinity (Mihalas 1978). Here 
the Sobolev optical thickness along the ray is defined by 



ts 



A L ' 



where L is the Sobolev length, 



L = v th /\v' 




(36) 



(37) 



and A is the width of the line, which in the case of Doppler 
broadening is 

^ _ / 3k B T M 

C V m atom 

With these definitions, equation ( |36| ) becomes 



(38) 



_ ^ijk 



(39) 



In the expanding Universe, the velocity gradient v' is given 
by the Hubble expansion rate H(z), and using the above 
definition for k. 



Ajj\j/_ [n^gj/gj) - nj] 
8nH(z) 



(40) 



To find the Sobolev escape probability for the ray, we 
average over the initial frequencies x, using the line profile 
function <fi(x) from equation (|29|), 



Pij 



dx<j)(x) exp 



~r s / 4>{x)dx = / d( exp(-rsC), 

J-oo J JO 



Prj 



1 - CXp(-Ts) 



(41) 



Note that this expression is independent of the line profile 
shape 4>(x). The escape probability p^ is defined as a fre- 
quency average at a single point. Finally, we must average 
over angle, but in the case of the isotropically expanding 
Universe, the angle-averaged Sobolev escape probability 
takes the same form as pij above. For further details on 
the Sobolev escape probability see Rybicki (1984) or Mi- 
halas (1978), §14.2. 

How does the Sobolev optical depth relate to the usual 
meaning of optical depth? The optical depth for a specific 
line at a specific redshift point is equivalent to the Sobolev 
optical depth, 



dr(uij) = -Ts(f>(x) — . 



(42) 



If no other line or continua photons are redshifted into that 
frequency range before or after the redshift point, then 
the optical depth at a given frequency today (i.e. summed 
over all redshift points) will be equivalent to the Sobolev 
optical depth at that past point. Generally, behaviors in 
frequency and space are interchangeable in a medium with 
a velocity gradient. 

In order to derive an expression for the bound-bound 
rate equations, we must consider the mean radiation field 
J in the line. For the case of spectral distortions J does 



Seager, Sasselov & Scott 



7 



not equal the Planck function at the line frequency. We 
use the core saturation method (Rybicki 1984) to get J 
using pij ; from this we get the net rate of deexcitations in 
that transition (j — > i), given in equation (p5|). In general, 
only the Lyman lines of H and He II, and the He I n 1 p-l 1 s 
lines have pij < 1. With this solution we have accom- 
plished two things: described the rcdshifting of photons 
through the resonance lines, and found a simple solution 
to the radiative transfer problem for all bound-bound lines. 

Peculiar velocities during the recombination era may 
cause line broadening of the same order of magnitude 
as thermal broadening over certain scales (A. Loeb, pri- 
vate communication). Because we use J(v) — B(i/,T), 
the radiation field is essentially constant over the width 
of the line and so the line shape is not important. Sim- 
ilarly, the peculiar velocities will not affect the Sobolev 
escape probability because it is independent of line shape 
(equation (|d|)). If peculiar velocities were angle depen- 
dent there would be an effect on the escape probability 
which is an angle-averaged function. Only for computing 
spectral distortions to the CMB, where the line shape is 
important, should line broadening from peculiar velocities 
be included. 

2.3.4. Chemistry 

Hydrogen molecular chemistry has been included be- 
cause it may affect the residual electron densities at low 
redshift (z < 200). During the recombination epoch, the 
H2 formation reactions include the H~ processes 



H + e~ • 
H" +H • 
and the H2 + processes 

H + H+ < 
H 2 + + H 

together with 



and 



H-+7, (43) 
H 2 +e", (44) 

H 2 ++7, (45) 
■H 2 +H+, (46) 

H + H 2 < — > H + H + H, (47) 

H 2 +e~< — >H + H + e". (48) 

(see Lepp & Shull 1984). The direct three-body process 
for H 2 formation is significant only at much higher densi- 
ties. The most recent rate coefficients and cross sections 
are listed with their references in Appendix A (see also 
Cen 1992; Puy et al. 1993; Tegmark et al. 1997; Abel et 
al. 1997; Galli & Palla 1998). 

We have not included the molecular chemistry of Li, 
He, or D. In general, molecular chemistry only becomes 
important at values of z < 200. Recent detailed analyses 
of H, D, and He chemistry (Stancil et al. 1996a) and of 
Li and H chemistry (Stancil et al. 1996b; Stancil & Dal- 
garno 1997) presented improved relative abundances of all 
atomic, ionic and molecular species. The values are cer- 
tainly too small to have any significant effect on the CMB 
power spectrum. 

There are two reasons for this. The visibility function 
is e~ T dT /dz, where r is the optical depth. The main com- 
ponent of the optical depth is Thomson scattering by free 



electrons. Because the populations of Li, LiH, HeH + , HD, 
and the other species are so small relative to the electron 
density, they do not affect the contributions from Thom- 
son scattering. Secondly, these atomic species and their 
molecules themselves make no contribution to the visibil- 
ity function because they have no strong opacities. HD 
has only a weak dipole moment. And while LiH has a 
very strong dipole moment, its opacity during this epoch 
is expected to be negligible because of its tiny (< 10~ 18 ) 
fractional abundance (Stancil, et al. 1996b). Similarly, the 
fractional abundance of H 2 D + (< 10 -22 ) is too small to 
have an effect on the CMB spectrum (Stancil et al. 1996a). 
For more details, see Palla et al. (1995) and Galli & Palla 
(1998). An interesting additional point is that because of 
the smaller energy gap between n=2 and the continuum in 
Li I, it actually recombincs at a slightly lower redshift than 
hydrogen does (see e.g. Galli & Palla 1998). Of course this 
has no significant cosmological effects. 

We have also excluded atomic D from the calculation. 
D, like H, has an atomic opacity much lower than Thomson 
scattering for the recombination era conditions. D par- 
allels H in its reactions with electrons and protons, and 
recombines in the same way and at the same time as H 
(see Stancil et al. 1996a). Although the abundance of D is 
small ([D/H] ~ 10~ 5 ), its Lyman photons are still trapped 
because they are shared with hydrogen. This is seen, for 
example, by the ratio of the isotopic shift of D Ly a to the 
width of the H Lya line, on the order of 10~ 2 . There- 
fore by excluding D we expect no change in the ionization 
fraction, and hence none in the visibility function. 

While the non-hydrogen chemistry is still extremely im- 
portant for cooling and triggering the collapse of primor- 
dial gas clouds, it is not relevant for CMB power spectrum 
observations at the level measurable by MAP and Planck. 

2.4. Expansion of the Universe 

The differential equations in time for the number densi- 
ties and matter temperature must be converted to differ- 
ential equations in redshift by multiplying by a factor of 
dz/dt . The redshift z is related to time by the expression 



^ = -(l + z)H(z), 



with scale factor 



a(t) = 



1 



1 + z 

Here H(z) = a/a is the Hubble factor 



(49) 



(50) 



H(z) 



2 _ tt2 
H 



fto(l+z) 4 /(l + ^cq) 



+ n (i + zf + njc(i + zf + n A , (51) 



where f2o is the density contribution, Qk is the curvature 
contribution and Qa is the contribution associated with 
the cosmological constant, with fio + + + = 1, 
and f^R = f2o/(l + z cq ). Here z cq is the redshift of matter- 
radiation equality, 



1 + z cq = O 



3(ctf ) 2 



87rG(l + /„)[/ : 



(52) 



with /„ the neutrino contribution to the energy density in 
relativistic species (/„ ~ 0.68 for three massless neutrino 
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types), G the gravitational constant, U the photon energy 
density, and H the Hubble constant today, which will be 
written as 100 h kms _:L Mpc . Since we are interested in 
rcdshifts z ~ z cq , it is crucial to include the radiation con- 
tribution explicitly (see HSSW). 

2.5. Matter Temperature 

The important processes that are considered in follow- 
ing the matter temperature are Compton cooling, adia- 
batic cooling, and Bremsstrahlung cooling. Less impor- 
tant but also included are photoionization heating, pho- 
torecombination cooling, radiative and collisional line cool- 
ing, collisional ionization cooling, and collisional recombi- 
nation cooling. Note that throughout the relevant time 
period, collisions and Coulomb scattering hold all the mat- 
ter species at very nearly the same temperature. 'Matter' 
here means protons (and other nuclei), plus electrons, plus 
neutral atoms; dark matter is assumed to be decoupled. 

Compton cooling is a major source of energy transfer 
between electrons and photons. It is described by the rate 
of transfer of energy per unit volume between photons and 
free electrons when the electrons are near thermal equilib- 
rium with the photons: 



dE t 



e,7 



4(JTUn c kB 



(It 



m c c 



(Tr - T M ), 



M 



dt 



8 cxt/rie 
3 rriecntot 



(Tr - T, 



(53) 



(54) 



(Weymann 1965), where E e> ~ is the electron energy den- 
sity, ks, m e and c are constants as before, ctt is the Thom- 
son scattering cross section, Tr is the radiation tempera- 
ture and Tm is the electron or matter temperature. To 
get from equation (53) to equation (^3) we use the energy 
of all particles; collisions among all particles keep them at 
the same temperature. Here ntot represents the total num- 
ber density of par ticles, which includes all of the species 
mentioned in §EJJ, while U represents the radiation en- 
ergy density (integrated over all frequencies) in units of 
ergs cm" 



-3. 



U = 



u(u)dv, 



(55) 



where u(v,t) = AirJ(i>,t)/c. In thermal equilibrium the 
radiation field has a frequency distribution given by the 
Planck function, J(v,t) — B(y,Tg), and thus, in thermal 
equilibrium the energy density is 



47T 

u(u,t) = — B V (T K ), 
c 



(56) 



and the total energy density U is given by Stefan's law 



U = 



The spectrum of the CMB remains close to blackbody be- 
cause the heat capacity of the radiation is very much larger 
than that of the matter (Peebles 1993), i.e. there are vastly 
more photons than baryons. 

Adiabatic cooling due to the expansion of the Universe 
is described by 



dT, 



M 



dt 



= -2H(t)Ti 



(58) 



since 7 = 5 /3 f° r an ideal gas implies Tm oc (1 + z) 2 . 
The following cooling and heating processes are often rep- 
resented by approximate expressions. We used the exact 
forms, with the exception of Bremsstrahlung cooling, and 
the negligible collisional cooling. 
Bremsstrahlung, or free-free cooling: 



A 



2 5 7re 6 Z 2 /27rfc B T 



brcm 



3 3 / 2 h 



P m c c° 



1/2 

gsn e (n p +nB_. 



cir 



-An 



Helll), 

(59) 

where gg is the free- free Gaunt factor (Seaton 1960), n p is 
the number density of protons, nneii and ?tHeiii the num- 
ber density of singly and double ionized helium respec- 
tively, and other symbols are as previously described. 
Photoionization heating: 



n^i = ^ UiAn I 
i=i Jv o 



cti C {v) 
hpv 



B(v,T R )h F (v - v )dv. (60) 



Photorccombination cooling: 

N , s LTE 



A r 



cuiv) 



n c n c 



hv>v 



^^ + B{v,T K ) 



- hpv l ksTw hj>{v-v )dv.(<al) 



Line cooling: 

Aime = hpvo[rijRji — riiRij] 
Collisional ionization cooling: 

A c _i = hpisoCi C . 
Collisional recombination heating: 

A c -r = hpi>oC C i. 



(62) 



(63) 



(64) 



Here vq is the frequency at the ionization edge. We used 
approximations for collisional ionization and recombina- 
tion cooling because these collisional processes are essen- 
tially negligible during the recombination era. Ci C and 
C c i are the collisional ionization and recombination rates 
respectively, computed as in e.g. Mihalas (1978), §5.4. 
Thus, with 

g?Tm dt dTu 

(65) 



dz 



dz dt 



the total rate of change of matter temperature with respect 
to redshift becomes 



(1 + *) 



2Tm 



dz 3H(z)m c c n c + nn + %e 

+ A c _i + A c _ r + Ai ino ) 



2(Ab rem Hp_i -f- Ap__ r 



3fc B n to t^(^) 



(T M - T R ) 

(66) 



Here nuo is the total number density of helium, and the 
denominator n G + nH + «Ho {fhat from equation (p4|)) takes 
into account the fact that the energy is shared among 
all the available matter particles. All the terms except 
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adiabatic cooling in equation ( pq ) involve matter energy 
conversion into photons. In particular, Compton and 
Bremsstrahlung cooling are the most important, and they 
can be thought of as keeping Tm very close to Tr until 
their time scales become long compared with the Hubble 
time, and thereafter the matter cools as Tm oc (1 + z) 2 . 
Previous recombination calculations only included Comp- 
ton and adiabatic cooling, however the additional terms 
add improvements only at the 10 _3 % level in the ioniza- 
tion fraction. The reason for the negligible improvement is 
that it makes little difference which mechanism keeps Tm 
close to Tr early on, and adiabatic cooling still becomes 
important at the same time. 

2.6. Summary of Equations 

The system of equations to be simultaneously integrated 
in redshift are: 



(1 + *) 



drii(z) 
dz 



7pr{ [n c (z)n c (z) 



H(z) 



Pd - ni(z)P ir 



E J N =1 Ai? JJ }+3n J (z), (67) 



(1 + *) 



dTi 



M 



8<ttU(J u , z ) 



dz 3H (z)m c c n c + n B + fift 



M 



Tr 



2T 



M 



2 (A t 



■n„ 



line } 



3k B n tot H(z) 



and 



(1 + *) 



(68) 



dJ(v, z) 
dz 



= 3J(v, z) — 



H{z) 



j(v,z)-K(u, Z)J(U, z) 



(69) 

For J{v,z) = B{v,z) = B(y,T K ) (see SgJ), equa- 
tion ( |69| ) can be omitted because the expansion of the Uni- 
verse preserves the thermal spectrum of non-interacting 
radiation, and we can use the Sobolev escape probability 
method for the primary spectral line distortions. The sys- 
tem of coupled equations (67) that we use contains up to 
609 separate equations, 300 for H (one for each of a max- 
imum of 300 levels we considered), 200 for He I, 100 for 
He II, 1 for He III, 1 for electrons, 1 for protons, and 1 
for each of the 5 molecular or ionic H species. This sys- 
tem of equations, along with ([38]), is extremely stiff, that 
is the dependent variables are changing on very different 
time scales. We used the Bader-Deufihard semi-implicit 
numerical integration scheme, which is described in Press 
et al. (1992). To test the numerical integration we checked 
at each time step that the total charge and total number 
of particles are conserved to one part in 10 7 . 

3. RESULTS AND DISCUSSION 

By an 'effective 3-level' H atom we mean a hydrogen 
atom that includes the ground state, first excited state, 
and continuum. In an effective 3-level atom, the energy 
levels between n=2 and the continuum are accounted for 
by a recombination coefficient which includes recombina- 
tions to those levels. This should be distinguished from 
an actual 3-level atom, which would completely neglect 
all levels above n=2, and would be a hopeless approxima- 
tion. Good accuracy is obtained by considering an n-level 
atom, where n is large enough. In practice we find that 



a 300-level atom is more than adequate. We do not ex- 
plicitly include angular momentum states £, whose effect 
we expect to be negligible. In contrast to the effective 3- 
level H atom, the 300-level H atom has no recombination 
coefficient with 'extra' levels. The 'standard' recombina- 
tion calculation refers to the calculation with the effective 
3-level atom that is widely used today and first derived by 
Peebles (1968) and Zel'dovich et al. (1968), updated with 
the most recent parameters and recombination coefficient 
(HSSW). 

The primordial He abundance was taken to be Yp = 0.24 
by mass (Schramm & Turner 1998). The present-day CMB 
temperature To was taken to be 2.728 K, the central value 
determined by the FIRAS experiment (Fixsen et al. 1996). 

3.1. The 'Effective 3-level' Hydrogen Atom 

For comparison with the standard recombination calcu- 
lation that only includes hydrogen (see Peebles 1968, 1993; 
Scott 1988), we reduce our chemical reaction network to an 
effective 3-level atom, i.e. a two-level hydrogen atom plus 
continuum. The higher atomic energy levels are included 
by way of the recombination coefficient, which can effec- 
tively include recombination to hundreds of levels. The 
following reactions are included: 



H n =2/=2s + 7 

H n= i + 7 
H n= i + 2 7 



e" +H+ 

H n= 2,^=2p 
H n= 2J=2s- 



As described in Peebles (1993), we omit the recombina- 
tions and photoionizations to the ground state because any 
recombination directly to the ground state will emit a pho- 
ton with energy > 13.6 eV, where there are few blackbody 
photons, and this will immediately re-ionize a neighboring 
H atom. We include the two-photon rate from the 2s state 
with the rate A 2s -i s = 8.22458s" 1 (Goldman 1989). The 
most accurate total Case B recombination coefficient is by 
Hummer (1994) and is fitted by the function 



qb = 10 



-13 



at" 



1 + ct d 



cm 3 s \ 



(70) 



where a = 4.309, b = -0.6166, c = 0.6703, d = 0.5300 and 
t = Tm/10 4 K (Pequignot et al. 1991; see also Verner & 
Ferland 1996). 

Consideration of detailed balance in the effective 3-level 
atom leads to a single ordinary differential equation for 
the ionization fraction: 

dx c _ [x 2 c n H a B - (3 B (l - x c )e-' lpl,2 °/ kBTM ] 



dz H(z)(l + z) 

[l + KA 2s -i s nB:(l - x c )] 

[l + K A 2s _i s n H (l - % e ) + K/3 B n B (l - x c )] 



(71) 



(see e.g. Peebles 1968; extra terms included in Jones & 
Wyse 1985, for example, are negligible). Here x c is the 
residual ionization fraction, that is the number of elec- 
trons compared to the total number of hydrogen nu- 
clei (»h). Here the Case B recombination coefficient 
as = ob(Tm), the total photoionization rate [3 B = 
a B (2Tr m c k B TM/ h 2 ,) 3 / 2 exp(— -E^sAbTm) as described in 
§2.3.2, V2s is the frequency of the 2s level from the ground 
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state, and the redshifting rate K = A^/ (8irH(z)), where 
X a is the Ly a rest wavelength. Note that Tm is used in 
equation ( |7l|) and in /3b, because the temperature terms 
come from detailed balance derivations that use Boltz- 
mann and Saha equilibrium distributions, which are colli- 
sional descriptions. In the past, this equation has been 
solve d ( for x e (z)) simultaneously with a form of equa- 
tion ( pq ) containing only adiabatic and Compton cooling. 
We refer to the approach of equation (|7l|) as the 'standard 
calculation.' 

For the comparison test with the standard recombina- 
tion calculation, we also use an effective recombination 
coefficient, but three equations to describe the three reac- 
tions listed above. That is, we simplified equation (|67]) to 
three equations, one for the ground state populationfni), 
one for the first excited state population (112), and one for 
the electrons (for H recombination n e — n p ). 



(1+*) 



dn\{z) 



1 



dz 



H{z) 



[Ai? 2p _ 



AR 



2s-ls 



-3ni (72) 



(1 + *) 



dn 2 (z) 
dz 



1 



(n c (z)n p (z)a B - n 2s {z)(3 B 

3n 2 (73) 



H(z) 

— AR 2 p-is — Ai? 2s _i 



(1 + ^) 



dn e (z) 
dz 



1 



H(z) 



[(n 2s (z)p B ~ n c (z)n p (z)a B ] + 3n 

The remaining physical difference between our effective 3- 
level atom approach and that of the standard calculation 
is the treatment of the redshifting of H Ly a photons (in- 
cluded in the A_R 2 p-is terms). In our calculation the red- 
shifti ng is accounted for by the Sobolev escape probability 
(see §pTf). Following Peebles (1968, 1993), the standard 



calculation accounts for the redshifting by approximating 
the intensity distribution as a step, and in effect takes the 
ratio of the redshifting of the photons through the line 
to the expansion scale that produces the same amount of 
redshifting. It can be shown that Peebles' step method 
considered as an escape probability scales as l/rg, where 
rs is the Sobolev optical depth. For high Sobolev optical 
depth, which holds for H Lya during recombination for 
any cosmological model (see Fig. 7), the Sobolev escape 
probability also scales as 1/rs: 



lim pij = lim — (1 

rs»l r s »l rs 



-TS\ _ \ 



(75) 



Therefore the two approximations are equivalent for Ly a, 
although we would expect differences for lines with rs <; 1, 
where p — > 1. Because we treat recombination in the same 
way as Peebles, no individual treatment of other lines is 
permitted, and therefore there are no other differences be- 
tween the two calculations for this simple case. Note that 
with Peebles' step method to compute AR 2p ^i s , and the 
assumption that ri\ — nu — n p , the above equations will 
reduce to the single ODE equation (|7l|). 

The results from our effective 3-level recombination cal- 
culation are shown in Fig. 1, plotted along with values 
from a separate code as used in HSSW, which represents 
the standard recombination calculation updated with the 
most recent parameters. The resulting ionization fractions 



are equal, which shows that our new approach gives ex- 
actly the standard result when reduced to an effective 3- 
level atom. Two other results are plotted for comparison, 
namely values of x c taken from Peebles (1968) and Jones 
& Wyse (1985). Their differences can be largely accounted 
for by the use of an inaccurate recombination coefficient 
with a B (T M ) oc T M 1/2 . 
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Fig. 1. — Comparison of effective three-level hydrogen recombi- 
nation for the parameters Otot = 1.0, f2e = 1-0, h = 1.0. Note that 
the Jones & Wyse (1985) and Peebles (1968) curves overlap as does 
our curve with the HSSW one. 

As an aside, we note the behavior for z <J 50 in our 
curve and the HSSW one. This is caused by inaccuracy 
in the recombination coefficient for very low temperatures. 
The down-turn is entirely artificial and could be removed 
by using an expression for a B (T) which is more physical 
at small temperatures. The results of our detailed calcu- 
lations are not believable at these redshifts either, since 
accurate modeling becomes increasingly difficult due to 
numerical precision as T approaches zero. But in any case 
the optical depth back to such redshifts is negligible, and 
the real Universe is reionized at a similar epoch (between 
z = 5 and 50 certainly). 

3.2. Multi-level Hydrogen Atom 

The purpose of a multi-level hydrogen atom is to im- 
prove the recombination calculation, by following the pop- 
ulation of each atomic energy level with redshift and by in- 
cluding all bound-bound and bound-free transitions. This 
includes recombination to, and photoionization from, all 
levels directly as a function of time, in place of a parame- 
terized recombination and photoionization coefficient. The 
individual treatment of all levels in a coupled manner al- 
lows for the development of departures from equilibrium 
among the states with time, and feedback on the rate of 
recombination. Since the accuracy of the recombination 
coefficient is probably the single most important effect in 
obtaining accurate power spectra (HSSW), it makes sense 
to follow the level populations as accurately as possible. 

In the multi-level H atom recombination calculation, we 
do not consider individual £ states (with the exception of 
2s and 2p), but assume the £ sublevels have populations 
proportional to (2£ + 1). The £ sublevels only deviate 
from this distribution in extreme non-equilibrium condi- 
tions (such as planetary nebulae). In their H recombi- 
nation calculation, Dell' Antonio & Rybicki (1993) looked 
for such £ level deviations for n < 10 and found none. For 
n > 10, the £ states are even less likely to differ from an 
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equilibrium distribution, because the energy gaps between 
the I sublevels are increasingly smaller as n increases. 




2000 



Fig. 2. — Multi-level hydrogen recombination for the standard 
CDM parameters Qtot = 1.0, S^b = 0.05, h = 0.5 (top), and for 
Otot = 1.0, n B = l.Q, h = 1.0 (bottom), both with Y P = 0.24, T = 
2.728 K. The 'effective 3-level' calculation is essentially the same as 
in HSSW, and uses a recombination coefficient which attempts to 
account for the net effect of all relevant levels. We find that we re- 
quire a model which considers close to 300 levels for full accuracy. 
Note also that although we plot all the way to z = 0, we know that 
the Universe becomes reionized at z > 5, and that our calculations 
(due to numerical precision at low T) become unreliable for z < 50 
in the upper two curves, and for z < 20 in the other curves. 



3.2.1. Results From a Multi-level H Atom 

Fig. 2 shows the ionization fraction x e from recombina- 
tion of a 2-, 10-, 50-, 100- and 300-level H atom, compared 
with the standard effective 3-level results. The x e con- 
verges for the highest n-level atom calculations. The effec- 
tive 3-level atom actually includes about 800 energy levels 
via the recombination coefficient (e.g. Hummer 1994). 
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Fig. 3. — Energy separations between various hydrogen atomic 
levels and the continuum. The solid curve shows the energy width of 
the same energy levels due to thermal broadening. Thermal broaden- 
ing was calculated using v(2k^Ty[/mnc 2 ) 1 /' 2 . The frequency of the 
highest atomic energy level (n=300, with an energy from the ground 
state of 109676.547 cm -1 ; the continuum energy level is 109677.766 
cm -1 ) was used for v, but a thermal broadening value for any atomic 
energy level would overlap on this graph. 

Fig. 2 shows that the more levels that are included in 
the hydrogen atom, the lower the residual x e . The simple 
explanation is that the probability for electron capture 
increases with more energy levels per atom. Once cap- 
tured, the electron can cascade downwards before being 
reionized. Together this means adding more higher energy 
levels per atom increases the rate of recombination. Even- 



tually x c converges as the atom becomes complete in terms 
of electron energy levels, i.e. when there is no gap between 
the highest energy level and the continuum (see Fig. 3). 
Ultimately the uppermost levels will have gaps to the con- 
tinuum which are smaller than the thermal broadening of 
those levels, and so energy levels higher than about n=300 
do not need to be considered, except perhaps at the very 
lowest redshifts. 

For other reasons entirely, our complete (300-level) H 
atom recombination calculation gives an x c lower than 
that of the effective 3-level atom calculation. The faster 
production of hydrogen atoms is due to non-equilibrium 
processes in the excited states of H, made obvious by our 
new, level-by-level tre atmen t of recombination. The de- 
tails are described in §3.2.2 below. 



3.2.2. 



Faster H Recombination in our Level-by-level 
Recombination Calculation 



The lower x c in our calculation compared to the stan- 
dard calculation is caused by the strong but cool radiation 
field. Specifically, both a faster downward cascade rate and 
a lower total photoionization rate contribute to a faster net 
recombination rate. 




Radiative rates n = 70 to n = 65 
Recombination rate to n = 70 
Radiative rates n = 70 to n = 60 
Radiative rates n = 70 to n = 20 
Radiative rates n = 70 to n = 2 



1200 1400 
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Fig. 4. — How the bound-bound rates of large energy separation 
(e.g. n=70 to n=20) go out of equilibrium at low T, illustrated us- 
ing an upper level of n=70 for definiteness. The case of equilibrium 
corresponds to net bound-bound rates of zero. Each upward and 
downward bound-bound rate for a given transition is represented by 
the same curve; non-equilibrium occurs where a single curve sepa- 
rates into two as redshift decreases. The rates shown are for the 
sCDM model. 

By following the population of each atomic energy level 
with redshift, we relax the assumption used in the stan- 
dard calculation that the excited states are in equilibrium. 
In addition, we calculate all bound-bound rates which con- 
trol equilibrium among the bound states. In the standard 
calculation, equilibrium among the excited states n > 2 
is assumed, meaning that the net bound-bound rates are 
zero. Fig. 4 shows that the net bound-bound rates are 
actually different from zero at z <; 1000. The reason for 
this is that at low temperatures, the strong but cool radi- 
ation field means that high energy transitions are rare due 
to few high energy photons. More specifically, photoex- 
citation and stimulated photodeexcitation for high energy 
transitions become rare (e.g. 70-10, 50-4 etc.). In this 
case spontaneous deexcitation dominates, causing a faster 
downward cascade to the n=2 state. In addition, the faster 
downward cascade rate is faster than the photoionization 
rate from the upper state, and one might view this as ra- 
diative decay stealing some of the depopulation 'flux' from 
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photoionization. Both the faster downward cascade and 
the lower photoionization rate contribute to the faster net 
recombination rate. 

The cool radiation field is strong, so photoexcitations 
and photodeexcitations are rapid among nearby energy 
levels (e.g. 70-65, etc.; see Fig. 4). What we see in Fig. 4 
is that with time, after z <; 1000, the n=70 energy level 
becomes progressively decoupled from the distant lower 
energy levels (n=2,3,...20,...), but remains tightly coupled 
to its nearby 'neighbors' (n=60,65, etc.). This explains 
the departures from an equilibrium Boltzmann distribu- 
tion (in the excited states) as seen in the shape of the 
curves in Fig. 5. 




principal quantum number n 




principal quantum number n 

Fig. 5. — Ratio of the actual number densities of the excited states 

(rii) and number densities for a Boltzmann distribution of excited 

states (n*) at different redshi fts for recombination within the stan- 

C — ^1 

dard CDM model. See section 3.2.2 for details. We give two different 
plots with linear and logarithmic y-axes and different redshift val- 
ues, to show a wider range of out-of-equilibrium conditions. The 
ratio approaches a constant for high n at a given z , because the high 
energy (Rydberg) states are all very close in energy, and thus have 
similar behavior (i.e. remain coupled to the radiation field and each 
other) . 

Fig. 5 illustrates the non-equilibrium of the excited 
states by showing the ratio of populations of the excited 
states compared to a Boltzmann equilibrium distribution 
with respect to n=2. We find that the upper levels of 
the hydrogen atom are not in thermal equilibrium with 
the radiation, i.e. the excited levels are not populated ac- 
cording to a Boltzmann distribution. The excited states 
are in fact overpopulated relative to a Boltzmann distri- 
bution. This is not a surprise for the population of the 
n=2 state, which is strongly overpopulated compared to 
the n=l ground state, and so should be all n > 2 states 
because all Lyman lines remain optically thick during re- 



combination. What is surprising is that all excited states 
develop a further overpopulation with respect to n=2 and 
each other. Note that this is not a population inversion. 
The recombination rate to a given high level is faster than 
the downward cascade rate, and this causes a 'bottleneck' 
creating the overpopulation. Fig. 5 shows that all states 
are in equilibrium at high redshifts, with the highest states 
going out of equilibrium first, followed by lower and lower 
states as the redshift decreases. The factor by which the 
excited states are over-populated approaches a constant 
at high n for a given redshift, with this factor increasing 
as z decreases. The ratio is constant because the high en- 
ergy level Rydberg states have very similar energy levels to 
each other, with a relatively large energy separation from 
the n=2 state (i.e. the exponential term in equation (21) 
dominates over the gi ratios, and the exponential term is 
similar for all of the Rydberg states.) 

Fig. 5 also shows an enormous ratio at low redshift 
(z < 500) for number densities of the actual excited states 
to the number densities of a Boltzmann distribution of 
excited states, on the order of 10 6 . At such a low red- 
shift, there are almost no electrons in the excited states 
(~ 10~ 20 cm -3 ), and so unlike at higher redshifts, the ra- 
tio is only an illustration of the strong departure from an 
equilibrium distribution; the actual populations are very 
low in any case. 

In comparison with the standard equilibrium capture- 
cascade calculation for cub , the unusual situation described 
above (caused by the strong but cool radiation field) leads 
to higher effective recombination rates for the majority of 
excited states without increasing photoionization propor- 
tionally. This results in a higher net rate of production of 
neutral hydrogen atoms, i.e. a lower x c . 

3.2.3. Accurate Recombination vs Recombination 
Coefficient 

To demonstrate why the non-equilibrium in the excited 
states of H affects the recombination rate, we must con- 
sider the difference in our new treatment of recombina- 
tion compared to the standard treatment. An important 
new benefit of our level- by- level calculation lies in replac- 
ing the recombination coefficient with a direct calculation 
of recombination to and photoionization from each level 
at each redshift step. In other words, we calculate the 
recombination rate and the photoionization rate using in- 
dividual level populations and parameters of the excited 
states i, 
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For the standard recombination calculation, the recombi- 
nation rate is 

N 
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and the photoionization rate is 

= n 2s a B (T M )e~ E ^/ k ^ (2^m c fc B T M ) 3/2 /h%. (79) 

In this last equation, the excited state populati ons a re hid- 
den by the Boltzmann relation with ri2s (see § [2.3.2 ) . The 
important point here is that our method allows redistribu- 
tion of the H level populations over all 300 levels at each 
redshift step, which feeds back on the recombination pro- 
cess via equation (|77j), and leads to the lower x c shown in 
Fig. 2. This redistribution of the level populations is not 
possible in the standard calculation's equation ( |79| ) which 
only considers the populations n c , m, and ri2 S , and consid- 
ers the excited level populations n > 2s to be proportional 
to n2 S in an equilibrium distribution. 

A small improvement in our new recombination treat- 
ment over the standard treatment is in our distinguishing 
the various temperature dependencies of recombination. 
Photoionization and stimulated recombination are radia- 
tive, so they should depend on Tr. Spontaneous recombi- 
nation is collisional and depends on Tm (see §2.3.1). In the 
standard calculation, the radiative nature of recombina- 
tion and photoionization is overlooked because both the re- 
combination coefficient and photoionization coefficient are 
a function Tm only (equations (J7q) and (|79j)). Although 
adiabatic cooling (equation (^)jdoes not aominate until 
quite low redshifts (z <J 100), it still contributes partially 
to matter cooling throughout recombination. The result- 
ing difference between Tm and Tr in the net recombination 
rate affects x c at the few percent level at z ^ 300 for the 
popular cosmologies, and has an even larger effect for high 
f^B models. 

3.2.4. Collisions 

The standard recombination calculation omits colli- 
sional excitation and ionization because at the relevant 
temperatures and densities they are negligible for a two- 
level hydrogen atom (Matsuda et al. 1971). We have found 
that the collisional processes are also not important for the 
higher levels, even though those electrons are bound with 
little energy. In high Qb models collisional ionization and 
collisional recombination rates for the highest energy levels 
are of the same order of magnitude as the photoionization 
and recombination rates, though not greater than them 
(see Fig. 6). 




Q n =1.0 n„=i.o H„=1.0 



Fig. 6. — Comparison of ionization rates. The upper curves 
are photoionization rates, the lower curves are collisional ionization 
rates. The left panel shows a model with the standard CDM param- 
eters, while the right panel shows an extreme baryon cosmology. 



3.2.5. Departures From Case B 

Case B recombination excludes recombination to the 
ground state and considers the Lyman lines to be opti- 
cally thick (i.e. photons associated with all permitted ra- 
diative transitions to n=l are assumed to be instantly re- 
absorbed). An implied assumption necessary to compute 
the photoionization rate is that the excited states (n > 2) 
are in equilibrium with the radiation. Unlike the standard 
recombination calculation, our method allows departures 
from Case B because the Lyman lines are treated by the 
Sobolev escape probability method which is valid for any 
optical thickness, and our method allow s departures from 
equilibrium of the excited states (§ 3.2.1 ). We find that the 
excited states depart from equilibrium at redshifts <J 800, 
so Case B does not hold then. However, our calculations 
show that for hydrogen all Lyman lines are indeed opti- 
cally thick during all of hydrogen recombination, so Case 
B holds for H recombination above redshifts ~ 800. 
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Fig. 7. — Sobolev optical depth (rg) for the H Lyman lines for 
two different cosmological models. From upper to lower the curves 
represent the optical depth in the Lyman transitions from n=2 
(i.e. Lyo), n=10, n=50, n=100, and n=270. The curves show that 
all the Lyman transitions are optically thick during H recombination 
(z < 2000) but some are optically thin (rg < 1) earlier, during He 

recombination (z > 2000). 

Fig. 7 shows that the Lyman lines are not optically thick 
at earlier times, e.g. during helium recombination, where 
we find some optically thin H Lyman lines. The Sobolev 
escape probability treats this consistently, which is neces- 
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sary because we evolve H, He I, and He II simultaneously, 
which is necessary because we evolve H, He I, and He II 
simultaneously. 

3.2.6. Other Recent Studies 

The previous study that was closest in approach to our 
own was that of Dell' Antonio & Rybicki (1993), who calcu- 
lated recombination for a ten-level hydrogen atom in order 
to estimate the spectral distortions to the CMB blackbody 
radiation spectrum. Ten levels arc insufficient to calculate 
recombination accurately, because the higher energy levels 
of the atom are completely ignored (see Fig. 3). However 
the accuracy of the ionization fraction (x G ) was sufficient 
to determine the magnitude of the spectral distortions. 
Their recombination model treated individual levels, but 
used a recombination coefficient to each level of the form 

1/9 

Tm • Because the form of the recombination coeffi- 
cient dominates the H recombination process, our models 
are not equivalent, and so there is little use in comparing 
the results. 

More recently, Boschan & Biltzinger (1998) derived a 
new parameterized recombination coefficient to solve the 
recombination equation of the standard calculation, and 
to generate spectral distortions in the CMB. Their cal- 
culation differs from ours in that their recombination co- 
efficient is pre-calculated. Hence it is not an interactive 
part of the calculation, and does not allow the advantages 
that our calculation does, mainly the feedback of the non- 
equilibrium in the excited states on the net recombination 
rate. While they include pressure broadening for a cutoff 
in the partition function, they neglect thermal broaden- 
ing. A more serious problem is their method of inclu- 
sion of stimulated recombination, as originally suggested 
by Sasaki & Takahara (1993), who included stimulated re- 
combination as positive recombination inst ead o f negative 
ionization. The physics (as described in §2.3.1) and our 
computational results are the same regardless of whether 
stimulated recombination is treated as positive recombi- 
nation or negative ionization. However, this may not be 
the case computationally for the standard calculation, if 
it is not treated with care. We defer a full discussion of 
these matters to a separate paper (Seager & Sasselov, in 
preparation) . 

We have also investigated how we can approximate our 
calculations, so that other researchers can obtain approx- 
imately accurate results without the need to follow 300 
levels in a hydrogen atom. Because the net effect of our 
new H calculation is a faster recombination (a lower freeze- 
out ionization fraction), our results can be reproduced by 
artificially speeding up recombination in the standard cal- 
culation. Further details are described in Seager, Sasselov, 
& Scott (1999). 

3.3. Helium 

We compute helium and hydrogen recombination simul- 
taneously. The recombination of He III into He II and 
He II into He I is calculated in much the same way as 
hydrogen, with recombination, photoionization, redshift- 
ing of the n 1 p-l 1 s lines (in H these are the Lyman lines), 
inclusion of the 2 1 s-l 1 s two-photon rates, collisional exci- 
tation, collisional deexcitation, collisiona l ioni z ation , and 
collisional recombination, as described in §2.3.1-1 



turn states up to the level n=20, above which only the 
principal quantum number energy levels and transitions 
are used. Fig. 8 (which shows the levels up to n=4 only) 
indicates how much more complicated the He I atom is 
compared with H or the hydrogenic He II. Our multi-level 
He II atom includes the first 4 angular momentum states 
up to the level n=4, above which only the principal quan- 
tum number energy levels and transitions are used. 
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Fig. 8. — Grotrian diagram for He I, showing the states with n < 4 
and the continuum. In practice our model atom explicitly contains 
the first 4 angular momentum states up to n=20, and 120 principal 
quantum number energy levels beyond. 

Photoionizations from any He I excited state are allowed 
only into the ground state of He II, because there are few 
photons energetic enough (> 40 eV) to do more than that. 
Two electron transitions in He I are negligible at recombi- 
nation era temperatures. 

Cosmological helium recombination was discussed ex- 
plicitly in Matsuda et al. (1969, 1971, hereafter MST), and 
Sato, Matsuda & Takeda (1971), and to a lesser extent in 
Lyubarsky & Sunyaev (1983), while several other papers 
give results, but no details (e.g. Lepp & Shull 1984; Fahr & 
Loch 1991; Galli & Palla 1998). The main improvement in 
our calculation over previous treatments of helium is that 
we use a multi-level He II atom, a multi-level He I atom 
with triplets and singlets treated correctly, and evolve the 
population of each energy level with redshift by includ- 
ing all bound-bound and bound-free transitions. This is 
not possible for the standard recombination calculation 
method (equation ([7l|)) extended to He I, using an effec- 
tive three-level He I atom with only a singlet ground state, 
singlet first excited state and continuum. 

3.3.1. Results From Hel Recombination 

Fig. 9, shows the ionization fraction x c through He II, 
He I and H recombination, plotted against the standard H 
calculation that includes He II and He I recombination via 
the Saha equation. For completeness wc give the helium 
Saha equations here: 
for He I ^ He II 



2.3.3 



_ The 

multi-level He I atom includes the first 4 angular momen- 
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and for He II <-> He III 



(a? e - 1 - /hcK = (27TTO G fc B r) / c - XHcII / fcB T (81) 
1 + 2/ Ho - /ipn H 

Here the %s are ionization potentials, nu is the total num- 
ber density of hydrogen, /n e is the total number fraction 
of helium to hydrogen /h c = jihc/«h = Yp /4(1 — Yp), and 
our definition of x c = n c /nn results in the complicated- 
looking left hand sides. The extra factor of 4 on the 
right hand side for He I <-* He II arises from the statistical 
weights factor. 
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Fig. 9. — Helium and hydrogen recombination for two cosmolog- 
ical models with Yp = 0.24 and T = 2.728 K. The first step from 
right to left is recombination of He III to He II, the second step is 
He II to He I, and the third step is H recombination. 
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Fig. 10. — Details of helium recombination for the standard CDM 
cosmology (top figure) and the high Ob cosmology (bottom figure). 
The dashed lines show our new results, and the dotted lines show 
the results assuming the 'standard calculation' (equivalent to Saha 
equilibrium). 

While our improved x e agr ees fa irly closely with Saha 
recombination for He II (see § 3.3.4 ), the difference in x c 
from Saha recombination during He I recombination is 
dramatic. Our new detailed treatment of He I shows He I 
recombination finishing just after the start of H recom- 
bination (see Fig. 9), i.e. significantly delayed compared 
with the Saha equilibrium case. This is different from the 
earlier calculations (e.g. MST), in which He I recombina- 
tion is finished well before H recombination begins. In 
this previous case, He I recombination still affected the 
CMB anisotropy power spectrum on small angular scales 



because the diffusion damping length grows continuously 
and is sensitive to the full thermal history (HSSW). In 
our new case, particularly for our low Qb models He I re- 
combination is still finishing at the very beginning of H 
recombination, which furthe r af fects the power spectrum 
at large angular scales (see §[3.7|). We show a 'blow-up' of 
the two helium recombination epochs in Fig. 10. 

3.3.2. Physics of He I Recombination 

The physics of He I recombination can be summarized 
as follows. There are three major aspects to it: (1) the 
He I has excited states which are able to retain charge; 
but (2) being very close to the continuum, the highly ex- 
cited states are easily photoionized by the radiation field 
at z ~ 3000; then (3) we have a standard hydrogenic-like 
Case B recombination, which is unaffected by neutral H 
removing He I 2 1 p-l 1 s (resonance line) photons. 
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Fig. 11. — : This figure shows why the Saha equilibrium recombi- 
nation rate (iJg a h a ) f° r H e I is n °t valid. Comparing the dotted and 
dashed lines, the photoexcitation rate (i.e. photoabsorption rate) for 
Hc 2 p— 1 s (i?Hc) is orders of magnitude greater than the photoion- 
ization rate for H (-Rh) from the same He I 2 1 p— l 1 s photon pool; 
there is no possibility for H to 'steal' the photons to speed up He I 
recombination. For Saha recombination to be valid, Rn > -Rhc, 
as well as > i?s a ha- For the sCDM model shown here, He I 
recombination begins around z = 3000. 

The He I atom has a metastable, i.e. very slow, set of 
states - the triplets (e.g. n 3 p-n 1 s). Therefore, overall the 
excited states of He I can naturally retain more charge 
than a simple hydrogenic system under Boltzmann equilib- 
rium. The situation would resemble what we found for H 
recombination with the enhanced populations of the higher 
states, and would lead to faster reduction of x c . However, 
the high excited states of He I are much more strongly 
'packed' towards the continuum compared to those of H; 
the energy difference between the 3p levels and the con- 
tinuum is 1.6 eV for He I versus 1.5 eV for H, compared to 
24.6 eV versus 13.6 eV for the ground state-continuum en- 
ergy difference. This is enough to depopulate the triplets 
(whose 'ground state' is n=2 s), given the much higher 
radiation temperature during He I recombination. Left on 
its own under these circumstances, He I would recombine 
much like the standard Case B effective 3-level H atom, 
i.e. slower than Saha recombination. There is one possi- 
ble obstacle - it is the existence of some neutral H, which 
could 'steal' He I resonance line photons, invalidate the ef- 
fective Case B and make it a Saha recombination instead. 
However, our detailed calculation shows that neutral H 
during He I recombination is not able to accomplish that, 
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and the process is not described by Saha equilibrium. 

Fig. 11 shows that Saha equilibrium recombination is in- 
valid for He I, by comparing the three possible destruction 
processes of the He I 2 1 p-l 1 s (in H this is Lya) photons: 
(1) cosmological redshift; (2) the 2 1 s-l 1 s two-photon rate; 
and (3) the photoionization rate of the ground state of H 
by the same He I 2 1 p-l 1 s photons. Fig. 11 clearly shows 
that process (3) is negligible (in contradiction to the dis- 
cussion in HSSW). To be doubly sure that absorption of 
these photons by hydrogen is negligible we explicitly in- 
cluded the relevant rate in our models and found no dis- 
cernible effects. In order for He I recombination to be 
approximated by Saha equilibrium, one of the 3 processes 
described above would have to be faster than or equal to 
the Saha equilibrium rate, which we do not find to be the 
case. 

The recombination of He I is slow for the same rea- 
sons that H recombination is, namely because of the op- 
tically thick n 1 p-l 1 s transitions which slow cascades to 
the ground state, and the exclusion of recombinations to 
the ground state. In other words He I follows a Case B 
recombination. Because the 'bottleneck' at n=2 controls 
recombination, it is not surprising that He I and H recom- 
bination occur at a similar redshift; the ionization energy 
of n=2 is similar in both. He I recombination is slower than 
H recombination because of its different atomic structure. 
The excited states of He I are more tightly packed, and 
the 2 1 p-l 1 s energy difference greater than that of H. The 
strong radiation field keeps the ratio of photoionization 
rate/downward cascade rate higher than in the H case, 
resulting in a slower recombination. 

We find the strong radiation field also causes the triplet 
states to be virtually unpopulated. The lack of electrons 
in triplet states is easily understood by considering the 
blackbody radiation spectrum. At He I recombination (z 
~ 3000), the blackbody radiation peak is around 2eV, so 
there are around 11 orders of magnitude more photons 
that can ionize the lowest triplet state 2 3 s (4.8 eV), than 
the singlet ground state (24.6 eV), since both are on the 
steeply decreasing Wien tail. It is interesting to note that 
in planetary nebulae where the young, hot ionizing star 
produces most of its energy in the UV, the opposite oc- 
curs: the He I atoms have few electrons in the singlet 
states; instead most of them are in the triplet states. 

There is one more possible method to speed up He I 
recombination, and that is collisional rates between the 
triplets and singlet states. If fast enough, the collisional 
rates would provide another channel to keep hold of cap- 
tured electrons - by pumping them into the triplet states 
faster than they can be reionized. The triplets are 3 times 
as populated as the singlets due to the statistical weight 
factors. By forcing the collisional rates to be greater than 
the recombination rates and the bound-bound radiative 
rates, we find an extremely fast He I recombination - ap- 
proximated by the Saha equilibrium. Essentially we force 
electrons from the singlets into the triplets faster than they 
can cascade downwards, and faster than they can be pho- 
toionized out of the triplets. In reality, the collisions are 
negligible, a few orders of magnitude less than the radia- 
tive rates. It is important to note that apart from colli- 
sions, the singlet and triplet states are only connected via 
the n 3 p-n 1 s transitions, which are orders of magnitude 
slower than the 2 3 s-l 1 s rate. We note here that MST 



stated that the collisional rates were high enough to cause 
equilibrium between the triplet and singlet states. One 
must be careful to compare all relevant rates, and we keep 
all of them in our code. We find the allowed radiative rates 
(e.g. photoexcitation and photodeexcitation) are greater 
than the collisional rates. Therefore the allowed radiative 
rates control the excited states' population distribution, 
not the collisional rates. In other words, electrons in the 
singlet states are jumping between bound singlet states 
faster than the collisional rates can send them into the 
triplet states. 

3.3.3. Effective 3-level calculation for He I 

We note here that MST used an effective 3-level He I sin- 
glet atom and calculated He I recombination in the same 
way as the standard H calculation (equation (|7l])) with the 
appropriate He I parameters. When we follow their treat- 
ment, we get essentially the same result as our multi- level 
He I calculation. We are not sure why MST obtained such 
a fast He I recombination. 

As with hydrogen, we have also investigated what is re- 
quired to achieve an accurate solution for helium, without 
modeling the full suite of atomic processes. We have found 
that the use of the 'effective 3-level' equations for helium 
(as described in MST), together with an appropriate re- 
combination coefficient for singlets only (equation (p2|)), 
results in a very accurate treatment of x c {z) duringthc 
time of helium recombination. In detail it is necessary 
to follow hydrogen and helium recombination simultane- 
ously, increasing the number of differential equations to 
solve. However, little accuracy is in fact lost by treating 
them independently - since recombination is governed by 
dramatic changes in time scales through Boltzmann fac- 
tors and the like, and is affected little by small changes 
in the number of free electrons at a given time. Further 
details are discussed in Seager, Sasselov, & Scott (1999). 

Although our model does not explicitly use a recombi- 
nation coefficient, it does allow us to calculate one easily. 
To aid other researchers it is worth presenting a fit for 
the singlet-only Case B recombination coefficient for He I 
(including recombinations to all states except the ground 
state) from the data in Hummer & Storey (1998). Hummer 
& Storey (1998) compute photoionization cross sections 
that are more accurate than the ones we use (Hofsaess 
1979), but are not publicly available. Following the func- 
tional forms used in the fits of Verner & Ferland (1996) we 
find 



QHc 




1+b 



m 3 s \ 



(82) 

with a = 10- 16,744 m 3 s- 1 , b = 0.711, T x = 10 5114 K, and 
T 2 fixed arbitrarily at 3 K. This fit is good to < 0.1% over 
the relevant temperature range (4,000-10,000 K), and still 
fairly accurate over a much wider range of temperatures. 

3.3.4. lie II Recombination 

He II recombination occurs too early to affect the power 
spectrum of CMB anisotropies. For completeness, we men- 
tion it briefly here. He II recombination is fast because 
of the very fast two photon rate. Fig. 14 shows that for 
most cosmologies the two-photon rate is faster than the 



Seager, Sasselov & Scott 



17 



net recombination rate, meaning that as fast as electrons 
are captured from the continuum they can cascade down 
to the ground state. Because of this, there is essentially 
no 'bottleneck' at the n=2 level. In high baryon models, 
He II recombination can be approximated using the Saha 
recombination. As shown in Fig. 9, He II recombination is 
slightly slower than the Saha recombination for low baryon 
models. 
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Fig. 12. — What controls H recombination? The net 2p-ls rate 
(dashed) compared to the 2s— Is two-photon rate (dotted) and the 
net recombination rate (solid) for 4 different cosmologies. Except 
for low Qb and low h models (e.g. the sCDM model), the 2s-ls rate 
dominates. The solid vertical line represents where 5% of the atoms 
have recombined. 



3.3.5. What Controls Recombination? 

H recombination is largely controlled by the 2s-ls two- 
photon rate, which except for low-baryon cases, is much 
faster than the H Ly a rate. The net recombination rate, 
net 2s-ls rate, and net Ly a rate are compared for dif- 
ferent cosmologies in Fig. 12. Figs. 13 and 14 show the 
same rate comparison for He I and He II. The three fig- 
ures all have the same scale on the x and y axes, for easy 
comparison. He I recombination is controlled by the 2 1 p- 
l 1 s rate rather than the 2 1 s-l 1 s rate as previously stated 
(e.g. MST). Fig. 13 (for He I) also illustrates the slow net 
recombination rate, which is the primary factor in the slow 
Case B He I recombination. Fig. 14 also illustrates that 
He II in the high fie and h models has a 2s-ls rate faster 
than the net recombination rate, meaning that there is 
no slowdown of recombination due to n= 2, and the Saha 
equilibrium approximation is valid. 

The rates change with cosmological model. Physically 
this is because all of the rates are very sensitive to the 
baryon density. The 2 1 p-l 1 s rates are further affected 
by the Hubble factor because the Sobolev approximation 
(equations ( ^l| ) and (|40|)), depends on the velocity gradi- 
ent. Whether most of the atoms in the Universe recom- 
bined via a 2p-ls or a 2s-ls two-photon transition de- 
pends on the precise values of the cosmological parameters. 
A confident answer to that question is still not known, 
given today's parameter uncertainties. 
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Fig. 13. — What controls He I recombination? The net 2 1 p-l 1 s 
rate (dashed) compared to 2 s— I s two-photon rate (dotted) and the 
net recombination rate (solid) for 4 different cosmologies. Except for 
high Qb and high h models, the 2 1 p— l J s rate dominates, in contrast 
to H. The solid vertical line represents where 5% of the atoms have 
recombined. 
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Fig. 14. — What controls He II recombination? The net 2p-ls 
rate (dashed) compared to the 2s-ls two-photon rate (dotted) and 
the net recombination rate (solid) for 4 different cosmologies. Except 
for low Qb an d l° w h models, the 2s-ls rate dominates during re- 
combination, and the 2p-ls at the start of recombination. The solid 
vertical line represents where 5% of the atoms have recombined. 



3.4. Atomic Data and Estimate of Uncertainties 

Our approach in this work has been to include all rel- 
evant degrees of freedom of the recombining matter in a 
consistent and coupled manner. This requires special at- 
tention to the quality of the atomic data used. The chal- 
lenge lies in building a consistent model for all energy lev- 
els and transitions, not just for the low-lying ones, which 
are often better known experimentally and theoretically. 

3.4.1. H and Hell 

Hydrogen (and the hydrogenic ion of helium) have 
exactly known rate coefficients for radiative processes 
from precise quantum-mechanical calculations (uncertain- 
ties below 1%). We use exact values for the bound-bound 
radiative transitions and for radiative recombination, as 
in e.g. Hummer (1994). For more details see Hummer 
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& Storey (1987), but also Brocklehurst (1970) and John- 
son (1972). In particular, the rate of radiative recombi- 
nation to level n of a hydrogenic ion can be evaluated 
from the photoionization (bound-free) cross section for 
level n, a nc (v ), wit h the standard assumption of detailed 
balance (see §2.3.1). For hydrogenic bound-free cross sec- 
tions, we follow in essence Seaton's work (Seaton 1959) 
with its asymptotic expansion for the Gaunt factor (see 
Brocklehurst 1970). Note that the weak dependence of 
the Gaunt factor on wavelength has a noticeable effect in 
our final recombination rate calculation. Given our ap- 
plication, we do not require the resolution of resonances, 
as achieved for a few transitions by the Opacity Project 
(TOPbase, Canto et al. 1993). Like Hummer (1994), we 
work with the n-levels assuming that the ^-sublevels have 
populations proportional to (2£+l). The resulting uncer- 
tainties for hydrogenic radiative rates at low temperatures 
(T < 10 5 K) certainly do not exceed the 1% level. 

Collisional rate coefficients cannot be calculated ex- 
actly. So, compared to the hydrogenic radiative rate coef- 
ficients, the situation for the bound-bound collisional rates 
and collisional ionization is poor, with errors typically 
about 6%, and as high as 20% in some cases (Percival & 
Richards 1978). A number of methods are used to evaluate 
electron-impact excitation cross sections of hydrogen-like 
ions (Fisher et al. 1997). These most recent values com- 
pare well to the older sources (Johnson 1972; Percival & 
Richards 1978). The helium ion, He II, is hydrogenic and 
was treated accordingly. We basically followed Hummer & 
Storey (1987) and Hummer (1994) in building the model 
atom. For our application, collisional processes are negli- 
gible, so that the large uncertainties that still persist for 
the collisional rates have no impact on our results. 

3.4.2. Eel 

Helium, in its neutral state, poses a challenge for build- 
ing a multi-level atomic model of high precision. Unlike 
atomic hydrogen, no exact solutions to the Schrodinger 
equation are available for helium. However, very high 
precision approximations are now available (Drake 1992, 
1994) which we have used. These approximations are es- 
sentially exact for all practical purposes. The largest rel- 
ativistic correction comes from singlet-triplet mixing be- 
tween states with the same n, L, and J, but is still small. 
Transition rates were calculated following the recent com- 
prehensive He I model built by Smits (1996) and some 
values in Theodosiou (1987). The source of our photoion- 
ization cross sections was TOPbase (Cunto et al. 1993) and 
Hofsaess (1979) for small n; above n=10 we used scaled 
hydrogenic values. New detailed calculations (Hummer 
& Storey 1998) show that the He I photoionization cross 
sections become strictly hydrogenic at about n > 20. The 
uncertainties in the He I radiative rates are at the 5% level 
and below. 

The situation with the collisional rates for He I is pre- 
dictably much worse than for He I radiative rates, with 
good R-matrix calculations existing only for n < 5 (Sawey 
& Berrington 1993). The collisional rates at large n are 
a crucial ingredient in determining the amount of singlet- 
triplet mixing, but fortunately collisions are not very im- 
portant for the low density conditions in the early Uni- 
verse, so the large uncertainty in these rates does not effect 
our calculation. The Born approximation, which assumes 



proportionality to the radiative transition rates, is used 
(see Smits 1996) to calculate the collisional cross sections 
for large n. 

For the 2 1 s-l 1 s two-photon rate for He I we used the 
value And = 51.3 s _1 (Drake et al. 1969) which differs 
from a previously used value (Dalgarno 1966) by ~ 10%. 
An uncertainty even of this magnitude would still make 
little difference in the final results. For the He II 2s-ls 
two-photon rate we used the value An e i = 526.5 s -1 (for 
hydrogenic ions this is essentially Z & times the value for 
H) from Lipeles et al. (1965). Dielectronic recombination 
for He I is not at all important during He I recombination. 
While dielectronic recombination dominates at tempera- 
tures above 6 x 10 4 K, for the range of temperatures rele- 
vant here it is at least 10 orders of magnitude below the 
radiative recombination rate (using the fit referred to in 
Abel et al. 1997). 

3.4.3. Combined Error From Atomic Data 

We have gathered together the uncertainties in the 
atomic data in order to estimate the resulting uncertainty 
in our derivation of x e . The atomic data with the dom- 
inant effect on our calculation are the set of bound-free 
cross sections for all H and He I levels - not so much any 
individual values, but the overall consistency of the sets 
(which are taken from different sources). The differences 
between our model atom and Hummer's (1994) reflect the 
uncertainty in the atomic data. To test the effect on our 
hydrogenic results, we compared the x c results of an effec- 
tive 3-level atom using Hummer's (1994) recombination 
coefficient with the results using a recombination coeffi- 
cient calculated with our own model H atom. We find 
maximum differences of 1% at z — 300, which corresponds 
to measurable effects on CMB anisotropies of much less 
than 1%. 

The error in x e arising from He I is more difficult to 
calculate. We estimate it to be considerably less than 1%, 
because the low level (n < 4) bound-bound and bound- 
free radiative rates dominate He I recombination, and as 
described above, those data are accurate. 

3.5. Secondary Distortions in the Radiation Field 

In our recombination calculation we follow 'secondary' 
distortions in the radiation field that could affect the re- 
combination process at a later time. The secondary distor- 
tions are caused by the primary distortions that are frozen 
into the radiation field. At a later time they are redshifted 
into interaction frequency with other atomic transitions. 
Explicitly, we follow: 

(1) H Ly a photons; 

(2) H 2s-ls photons. 

By the time of H recombination these photons have been 
redshifted into an energy range where they could photoion- 
ize H(n=2). In addition we follow: 

(3) He I 2 1 p-l 1 s; 

(4) He I 2 1 s-l 1 s. 

By the time of H recombination these photons have been 
redshifted into an energy range which could photoionizc 
H(n=l). And finally we also follow: 

(5) He II Ly a photons; 

(6) He II 2s-ls photons. 

By the time of H recombination these photons have been 
redshifted into an energy range which could photoionize 
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H(n=l). These He II photons bypass He I because the 
photons have not been redshifted into a suitable energy 
range for interaction. 

Here we only attempt to investigate the maximum ef- 
fects of secondary spectral distortions. To that end we do 
not include additional distortions which are smaller. For 
example, Lyman lines other than (1), (3), and (5), whose 
distortions are smaller than Ly a, will produce a compa- 
rably smaller feedback on photoionization. The He I sin- 
glet recombination photons could theoretically photoionize 
He I triplet states, but as previously discussed there are 
virtually no electrons in the triplet states, so this process is 
also negligible. Another possible effect is due to the similar 
energy levels of H and He II: Ai?Hcii = 4A£"h- For exam- 
ple, the transition from He II (n=4) to (n=2) produces the 
same frequency photons as the transition from H (n=2) 
to (n=l). These transitions are theoretically competing 
for photons, and this effect can be important for other 
astrophysical situations (e.g. planetary nebulae) where H 
and He II simultaneously exist. However, any such ef- 
fect is negligible for primeval recombination because dur- 
ing He II recombination the amount of neutral H is very 
small ([H/He II] < 10~ 8 ), and during H recombination, 
there is almost no He II ([He II/H]< 1CT 10 ). 

Because we are only investigating maximum effects we 
assume the photons were emitted at line center and are 
redshifted undisturbed until their interaction with H(n=l) 
or H(n=2) as described above. We also assume two pho- 
tons at half the energy for the 2s-ls transitions, compared 
to the 2p-ls transitions. The distorting photons emitted 
at a time z em are absorbed at a later time z, where 



Z Zcm^edge/^e 



(83) 



Here f e dge is the photoionization edge frequency where the 
photons are being absorbed, and v em is the photon's fre- 
quency at emission. The distortions are calculated as 



J(v,z) = hi/(z)cpij(z em ) x 



rii{z em )BijB{y era ,T R {z em )) V, (84) 



where: B(v em ,Tji(z em )) is the Planck function at the time 
of emission; Aji, Bji and Bij are the Einstein coefficients; 
Pij is the Sobolev escape probability for the line; and the 
other variables are as described previously. 

The distortions (1) and (2) were previously discussed by 
Rybicki & Dell' Antonio (1992). They pointed out that the 
effect from (1) should be small, because the Lya distor- 
tion must be redshifted by at least a factor of 3 to have 
any effect. This means that the Ly a photons produced at 
z <J 2500 will only affect the Balmer continuum at z <J 800 
when the recombination process (and any possibility of 
photoionization) is almost entirely over. 

We find that including the distortions (1) through (6) 
improves x c during H recombination at less than the 0.01% 
level. This difference is far too small to make a significant 
change in the power spectrum, and it is negligible com- 
pared to the major improvements in this paper, which are 
the level-by-level treatment of H, He I, and He II, allowing 



departures of the excited state populations from an equi- 
librium distribution, calculating recombination directly to 
each excited state, and the correct treatment of He I triplet 
and singlet states. However, the removal of these distort- 
ing photons by photoionization must be taken into account 
when calculating spectral distortions to the CMB black- 
body, which we plan to study in a later paper. 

3.6. Chemistry 



Including the detailed hydrogen chemistry (see §2.3.4) 



marginally affects the fractional abundances of protons 
and electrons at low z. However, the correction is of the 
order 10 _2 x c at z < 150. This change in the electron den- 
sity would change the Thomson scattering optical depth 
by ~ 10~ 5 , too little to make a difference in the CMB 
power spectrum. 




200 220 240 260 280 300 

redshift 

Fig. 15. — The effect of the improved treatment of recombination 
on H chemistry. Shown is the standard CDM model. Solid lines are 
values from the standard calculation, dashed and dotted lines from 
our improved results. 

On the other hand, as shown in Fig. 15, the different 
x e (z) that we find will lead to fractional changes of sim- 
ilar size in molecular abundances at low z, since H2 for 
example is formed via H~ which is affected by the resid- 
ual free electron density. The delay in He I recombination 
compared to previous studies causes a similar delay in for- 
mation of He molecules (P. Stancil, private communica- 
tion). However, with the exception of He J, no changes 
are greater than those caused by the residual free electron 
density at freeze-out. Since molecules can be important for 
the cooling of primordial gas clouds and the formation of 
the first objects in the Universe, the precise determination 
of molecular abundances is an important issue (e.g. Lepp 
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& Shull 1984, Tegmark et al. 1997, Abel et al. 1997, Galli 
& Palla 1998). However, the roughly 10-20% change in the 
abundance of some chemical species is probably less than 
other uncertainties in the reaction rates (A. Dalgarno, pri- 
vate communication). With this in mind, we suspect no 
drastic implications for theories of structure formation. 
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Fig. 16. — Differences in CMB power spectra arising from the 
improved treatment of hydrogen for: (a) the standard CDM param- 
eters Qtot = 10, Q B = 0.05, h = 0.5, Y P = 0.24, T = 2.728 K; 
and (b) an extreme baryon model with Qb = 1, h = 1.0 (for which 
there is relatively little effect). The fractional difference plotted is 
between our new hydrogen recombination calculation and the stan- 
dard hydrogen recombination calculation (e.g. in HSSW), with the 
sense of C° ew — C£ ld , and with the two calculations normalized to 
have the same amplitude for the initial conditions. The solid lines 
are for temperature, while the dashed lines are for the ('E'-modc) 
polarization power spectrum. 



3.7. Power Spectrum 

Even relatively small differences in the recombination 
history of the Universe can have potentially measurable 
effects on the CMB anisotropics. And so we might expect 
our two main changes (one in H and one in He) to be no- 



ticeable in the power spectrum. As a first example Fig. 16 
compares the difference in the anisotropy power spectrum 
derived from our new x c (z) to that derived from the stan- 
dard recombination x c (essentially identical to that de- 
scribed in HSSW) , for hydrogen recombination only. Here 
the Cgs are squares of the amplitudes in a spherical har- 
monic decomposition of anisotropies on the sky (the az- 
imuthal index m depends on the choice of axis, and so 
is irrelevant for an isotropic Universe). They represent 
the power and angular scale of the CMB anisotropies by 
describing the rms temperatures at fixed angular separa- 
tions averaged over the whole sky (see e.g. White, Scott, 
& Silk 1994). These Cgs depend on the ionization frac- 
tion x c through the precise shape of the thickness of the 
photon last scattering surface (i.e. the visibility function). 
Since the detailed shape of the power spectrum may al- 
low determination of fundamental cosmological parame- 
ters, the significance of the change in x c is evident. To 
determine the effect of the change in x c we have used the 
code cmbf ast written and made available by Seljak & Zal- 
darriaga (1996), with a slight modification to allow for the 
input of an arbitrary recombination history. 

The dominant physical affect arising from the new H 
calculation comes from the change in x c at low z. A 
process seldom mentioned in discussions of CMB aniso- 
tropy physics (which are otherwise quite comprehensive, 
e.g. Hu, Silk, & Sugiyama 1997) is that the \aw-z tail 
of the visibility function results in partial erasure of the 
anisotropies produced at z ~ 1000. The optical depth 
in Thomson scattering back to, for example, z — 800 
(r = cctt / n e (dt/dz) dz) can be several percent. This par- 
tial rescattering of the photons leads to partial erasure of 
the Cgs by an amount e _2r . Let us look at the standard 
CDM calculation first (Fig. 16(a)). Our change in the op- 
tical depth back to z ~ 800 (see Fig. 2) is around 1% less 
than that obtained using the standard calculation, and so 
we find that the anisotropies suffer less partial erasure by 
about 2%. There is no effect on angular scales larger than 
the horizon at the scattering epoch (here redshifts of sev- 
eral hundred), so that all multipoles are effected except 
for the lowest hundred or so is. Hence this effect is largely 
a change in the overall normalization of the power spec- 
trum, with some additional differences at low I which will 
be masked by the 'cosmic variance'. In addition there are 
smaller effects due to changes in the generation of anisotro- 
pies in the low-z tail, giving small changes in the acoustic 
peaks, which can be seen as wiggles in the figure. Since the 
partial erasing effect is essentially unchanged in the case 
of the Ob = h = 1.0 model, these otherwise sub-dominant 
effects are more obvious in Fig. 16(b). 

Differences in the power spectra are rather small in ab- 
solute terms, so Fig. 16 plots the relative difference. We 
have shown this for our two chosen models, one being 
standard Cold Dark Matter (a), which we will refer to 
as sCDM, and the other being an extreme baryon-only 
model (b). These models are meant to be representative 
only, and changes in cosmological parameters will result 
in curves which differ in detail. We describe how to cal- 
culate an approximately correct recombination history for 
arbitrary models in a separate paper (Seager, Sasselov, & 
Scott 1999). Since the main effect is similar to an overall 
amplitude change, we normalized our CMB power spectra 
to have the same large-scale matter power spectrum, which 
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is equivalent to normalizing to the same amplitude for the 
initial conditions. The amplitude of the effect of our new H 
calculation clearly depends on the cosmology. For the high 
fie and h case, the freeze-out value of x c is much smaller 
(around 10~ 5 ), and since the fractional change in x e is sim- 
ilarly ~ 10%, the absolute change in ionization fraction is 
much lower than for the sCDM model. The integrated op- 
tical depth is directly proportional to fifth Ax c , which is 
small, despite the increase in fie and h. Hence we see a 
much smaller increase from our hydrogen improvement in 
Fig. 16(b). The normalization change is rather difficult to 
see, since it is masked by relatively small changes around 
the power spectrum peaks, giving wiggles in the difference 
spectrum. 

The dashed lines in Fig. 16 show the effect on the power 
spectrum for CMB polarization. In standard models po- 
larization is typically at the level of a few percent of the 
anisotropy signal, and so will be difficult to measure in 
detail (see Hu & White 1997b for a discussion of CMB po- 
larization) . We show the results here to indicate that there 
are further observational consequences of our improved re- 
combination calculation (explicitly we have plotted the 'E' 
mode of polarization, see e.g. Seljak 1997). The effect of 
our improvements on the polarization can be understood 
similarly through the visibility function. Since the po- 
larization power spectrum tends to have sharper acous- 
tic peaks, the wiggles in the difference spectrum are more 
pronounced than for the temperature anisotropies. Note 
that the large relative differences at low I are actually 
very small in absolute terms, since the polarization signal 
is so small there. The polarization-temperature correla- 
tion power spectrum and the 'B' mode of polarization (for 
models with gravity waves) could also be plotted, but little 
extra insight is gained, and so we avoid this for the sake 
of clarity. 

The other major difference we find compared with pre- 
vious treatments is in the delayed recombination of He I. 
In Fig. 17 we show the effect of our new He I calculation, 
again as a fractional change in the CMB anisotropy power 
spectrum versus multipole I. The change in the recom- 
bination of He I affects the density of free electrons just 
before hydrogen recombination, which in turn affects the 
diffusion of the photons and baryons, and hence the damp- 
ing scale for the acoustic oscillations which give rise to the 
peaks in the power spectrum. The phases of the acoustic 
oscillations will also be affected somewhat, which shows 
up in the wiggles in the difference spectrum. For CDM- 
like models the main effect is the change in the damping 
scale, since we now think there are more free electrons at 
z <~ 1500-2000. The resulting change in the Cgs is essen- 
tially the same as assuming the wrong angular scale for 
the damping of the anisotropies (see Hu & White 1997a), 
which is the same physical effect that HSSW found in ar- 
guing for the need to include He I recombination at all 
for obtaining percent accuracy in the Cgs. The effect of 
this improved He I on the power spectrum will depend 
on the background cosmology through the baryon density 
(oc flBh 2 ) and the horizon size at last scattering through 
floh 2 - hence there is no simple fitting formula, and it is 
necessary to calculate the effect on the anisotropy damping 
tail for each cosmological model considered. 
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Fig. 17. — The effect of the improved treatment of helium on the 
CMB power spectra for: (a) the standard CDM model; and (b) the 
extreme baryonic model (for which there is essentially no change). 
The fractional difference plotted is between our new helium recom- 
bination calculation and the assumption that helium follows Saha 
equilibrium (as in HSSW), with the same 'effective 3-level' hydrogen 
recombination used in both cases. Again the sense is C^ ew — C£ ld , 
solid lines are temperature, and dashed lines are polarization. 

There are really two parts to the He I effect. Firstly the 
extra x c makes the tight coupling regime tighter, so that 
the photon mean free path is shorter, and the length scale 
for diffusion is smaller. Secondly, the effective damping 
scale comes from an average over the visibility function, so 
an increase in the high-z tail also leads to a smaller damp- 
ing scale. The CMB anisotropies can be thought of as a 
series of acoustic peaks multiplied by a roughly exponen- 
tial damping envelope, with the characteristic multipole of 
the cut-off being determined by the damping length scale. 
As a result of this smaller damping scale, the high i part 
of the power spectrum is less suppressed, and so we see an 
increase in Fig. 17(a) towards high I. For the CIb = h = 1 
model (Fig. 17(b)) we see only a very small effect at the 
highest is. This is easily understood by examining Fig. 10, 
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where we see that He I recombination is pushed back to 
higher rcdshifts than for the sCDM case, and also that H 
recombination happens earlier, which, together with the 
higher £1b and h, shifts the peak of the visibility func- 
tion to lower redshifts relative to the recombination curve. 
Hence the the high-z tail of the visibility function is much 
less affected in this case, and our improved He I calcu- 
lation has essentially negligible effect. However, for less 
extreme models we find that the He I effect is always at 
least marginally significant. 
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Fig. 18. — The total effect of our improvements on the CMB 
power spectra for: (a) the standard CDM model; and (b) an ex- 
treme baryonic model. These plots are essentially the sum of the 
separate effects of hydrogen and helium. 

Taking the two main effects together, for the sCDM 
model, we find that they have essentially the same sign, 
so that the total effect of our new calculation is more dra- 
matic (see Fig. 18(a)). Both effects lead to a slight in- 
crease in the anisotropies, particularly at small angular 
scales (high £) . Although the exact details will depend on 
the underlying cosmological model, we see the same gen- 
eral trend for most parameters we have looked at. The 



change can be higher than 5% for the smallest scales, al- 
though it should be remembered that the absolute ampli- 
tude of the anisotropies at such scales is actually quite low. 
Nevertheless, changes of this size are well above the level 
which is relevant for future determinations of the power 
spectrum. As a rough measure, the cosmic variance at 
I ~ 1000 is about 3%. Hence a 1% change over a range 
of say 1000 multipoles is something like a 10cr effect for a 
cosmic variance limited experiment. What we can see is 
that although the effects are far from astonishing, they are 
at a level which is potentially measurable. Hence our im- 
provements are significant in terms of using future CMB 
data-sets to infer the values of cosmological parameters; if 
not properly taken into account these subtle effects in the 
atomic physics of hydrogen and helium might introduce 
biases in the determination of fundamental parameters. 

3.8. Spectral Distortions to the CMB 

With the model described in this paper we plan to cal- 
culate spectral distortions to the blackbody radiation of 
the CMB today (see also Dubrovich 1975; Lyubarsky & 
Sunyaev 1983; ' Fahr & Loch 1991; Dell' Antonio & Ry- 
bicki 1993; Burdyuzha & Chekmezov 1994; Dubrovich & 
Stolyarov 1995, 1997; Boschan & Biltzinger 1998). The 
main emission from Ly a and the two- photon process will 
be in the far-infrared part of the spectrum. Transitions 
among the very high energy levels, which have very small 
energy separations, may produce spectral distortions in 
the radio. Although far weaker than distortions by the 
lower Lyman lines, they will be in a spectral region less 
contaminated by background sources. 

Detecting such distortions will not be easy, since they 
are generally swamped by Galactic infra-red or radio emis- 
sion and other foregrounds. Our new calculation does not 
yield any vast improvement in the prospects for detection. 
However, confirmation of the presence of these recombina- 
tion lines would be a definitive piece of supporting evidence 
for the whole Big Bang paradigm. Moreover, detailed mea- 
surement of the lines, if ever possible to carry out, would 
be a direct diagnostic of the recombination process. For 
these reasons we will present spectral results elsewhere. 

4. CONCLUSIONS 

One point we would like to stress is that our detailed cal- 
culation agrees very well with the results of the effective 
3-level atom. This underscores the tremendous achieve- 
ment of Peebles, Zel'dovich and colleagues in so fully un- 
derstanding cosmic recombination 30 years ago. However, 
the great goal of modern cosmology is to determine the 
cosmological parameters to an unprecedented level of pre- 
cision, and in order to do so it is now necessary to under- 
stand very basic things, like recombination, much more 
accurately. 

We have shown that improvements upon previous re- 
combination calculations result in a roughly 10% change 
in x c at low redshift for most cosmological models, plus 
a substantial delay in He I recombination, resulting in a 
few percent change in the CMB power spectrum at small 
angular scales. Specifically, the low redshift difference in 
x e is due to the H excited states' departure from an equi- 
librium distribution. This in turn comes from the level-by- 
level treatment of a 300-level H atom, which includes all 
bound-bound radiative rates, and which allows feedback 
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of the disequilibrium of the excited states on the recombi- 
nation process. The large improvement in x c during He I 
recombination comes from the correct treatment of the 
atomic levels, including triplet and singlet states. While 
it was already understood that He I recombination would 
affect the power spectrum at high multipoles (HSSW), our 
improved He I recombination affects even the start of H re- 
combination for traditional low 17b models. There is thus 
a substantially bigger change in the C^s, reaching to larger 
angular scales. 

Careful use of Tm rather than Xr can also have notice- 
able consequences, as to a lesser extent can the treatment 
of Lya redshifting using the Sobolev escape probability. 
Our other new contributions to the recombination calcula- 
tion produce negligible differences in x c . Collisional exci- 
tation and ionization for H, He I and for He II are of little 
importance. Inclusion of additional cooling and heating 
terms in the evolution of Tm also produce little change in 
x e . The largest spectral distortions do not feed back on the 
recombination process to a level greater than 0.01% in x c . 
Finally, the H chemistry occurs too low in redshift to make 
any noticeable difference in the CMB power spectrum. 

Although we have tried to be careful to consider every 



process we can think of, it is certainly possible that other 
subtle effects remain to be uncovered. We hope that we 
do not have to wait another 30 years for the next piece 
of substantial progress in understanding how the Universe 
became neutral. 

The program recfast, which performs approximate 
calculations of the recombination hist ory is available at 



http : //www. astro . ubc . ca/people/scott/recf ast .html 



(FORTRAN version) and at jhttp : //cf a-www . harvard . eduA 
^sasselov/rec/ (C version^ We would like to thank 
George Rybicki, Ian DeU'Antonio, Avi Loeb, and Han 
Uitenbroek for many useful conversations. Also David 
Hummer and Alex Dalgarno for discussions on the atomic 
physics, Alex Dalgarno and Phil Standi for discussions on 
the chemistry, Martin White, Wayne Hu and Uros Seljak 
for discussion on the cosmology, and Jim Peebles for dis- 
cussions on several aspects of this work. We thank the 
referee for a careful reading of the manuscript. Our study 
of effects on CMB anisotropics was made much easier 
through the availability of Matias Zaldarriaga and Uros 
Seljak's code cmbf ast. DS is supported by the Canadian 
Natural Sciences and Engineering Research Council. 
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APPENDIX 

REACTION RATES AND CROSS SECTIONS 
Where cross sections are listed instead of reaction rates, we calcul ate th e reaction rate through the integrals for photoion- 



ization (or photodissociation) and recombination as described in §2.3.1. 

[1] Dc Jong 1972; [2] Hirasawa 1969; [3] Karpas et al. 1979; [4] Shapiro & Kang 1987; [5] Dove & Mandy 1986; [6] Gredel 
& Dalgarno 1995; [7] Wishart 1979 and Broad & Reinhardt 1976. 



Reaction 




Rate Coefficient (cm 3 s L ) 


Reference 


H- + H — 


-+ H 2 + e- 


1.30 x 10" y 


[1] 


H 2 + e~ - 


H + H 


2.70 x 10- 8 T^ 3/2 exp(-43000/T M ) 


[2] 


H 2 + + H - 


-> H 2 + H+ 


6.40 x 10~ 10 


[3] 


H 2 + H+ - 


-> H 2 + + H 


2.4 x 10- 9 exp(-21200/T M ) 


W 


H + H 2 <— 


-> H + H + H 


1.0 x 10~ 10 exp(-52000/r M ) 


[5] 


H 2 + e" <- 


-> H + H + e" 


2.0 x 10~ 9 (TM/300) - 5 exp(-116300/r M ) 


[6] 



Reaction 



H + H+ 

-> H + e- 
H+ + e- 
+ He+ + e" 
He++ + e" 



Cross Section (cm 2 ) 



Reference 



H 2 ++ 7 

H~ + 7^ 
H + 7 «- 
He + 7 <- 
He+ + 7 



See expression in reference. 
See table in reference. 



See section 3.4 
See section 3.4 
See section 3.4 



[4] 
[7] 



